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INTRODUCTION 


During  the  past  four  years,  under  the  support  of  Office  of  Naval  Research,  a 
research  effort  has  been  underway  on  the  following  two  topics : 

(a)  Numerical  simulation  of  supersonic  shear  layer  mixing  phenomena. 

(b)  Development  of  efficient  methods  for  3-D  unsteady  incompressible  viscous  flow 
simulations. 

Substantial  progress  was  made  in  the  above  areas,  and  the  computer  codes 
developed  as  part  of  these  efforts  were  transfored  for  further  use  to  Virgiiua  Polytechnic 
Institute  (Dr.  Saad  Ragab)  and  to  David  Taylor  Research  and  Development  Center  (Dr.  Wei 
Tang). 

The  progress  made  has  been  documented  in  a  number  of  conference  proceedings, 
two  journal  articles  and  a  Ph.  D.  dissertation  ( in  progress).  The  Appendix  contains  all  the 
published  work  done  to  date. 

This  research  led  to  a  number  of  significant  new  findings.  Some  of  the  major 
accomplishments  are  as  follows : 

(a)  A  highly  accurate  method  for  simulating  compressible  mixing  layers  was  developed. 
This  method  is  suitable  for  direct  numerical  simulation  (DNS)  and  large  eddy  simulation 
(LES)  of  compressible  turbulence,  and  is  used  in  this  context  by  other  researchers. 

(b)  An  iterative  time  marching  method  for  3-D  incompressible  flows  was  developed.  This 
method  is  robust,  and  can  handle  internal  and  external  flows.  It  employs  a  multigrid 
iterative  strategy  for  satisfying  the  discretized  form  of  the  governing  equations  of  3-D 
viscous  flows  to  great  accuracy. 


It  is  hoped  that  these  efforts  will  serve  as  useful  stepping  stones  for  future  research 


in  these  areas. 
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ABSTRACT 

The  behavior  of  free  shear  layers  within 
ramjet  dump  commbustors  is  studied  through  the 
numerical  solution  of  unsteady  compressible 
Navier-Stokes  equations.  Three  configurations  are 
considered:  a)  a  short  combustor  with  an  open 
downstream  boundary,  b)  a  long  combustor  with  an 
open  downstream  boundary,  and  c)  a  short  combustor 
with  a  partially  blocked  downstream  boundary. 
Vorticity  contours  of  the  computed  flow  fields  in 
a11  the  three  cases  reveal  oscillations  of  the 
shear  layer,  roll  up  and  shedding  of  organized 
vortices.  A  Fourier  analysis  of  the  computed  flow 
fields  indicates  that  the  natural  acoustic  fre¬ 
quency  of  the  system,  and  the  natural  shear  layer 
instability  frequency  are  the  two  dominant  fre¬ 
quencies  of  the  flow  field.  It  is  also  observed 
that  the  boundary  conditions  play  a  crucial  role 
in  the  behavior  of  the  combustor  flow  field. 

INTRODUCTION 

The  flow  field  within  ramjet  combustors  is 
characterized  by  a  variety  of  phenomena  such  as 
thin  boundary  layers  along  the  walls  of  the  inlet 
and  combustor,  recirculating  flow  zones,  reacting 
flow  and  free  shear  layers.  The  free  shear  layer 
is  unstable  by  nature,  and  can  undergo  large 
spatial  and  temporal  oscillations  when  subjected 
to  disturbances.  In  many  cases,  the  shear  layer 
may  roll  up  and  form  vortices,  which  are  shed  at 
periodic  intervals  resulting  in  a  highly  unsteady, 
and  undesirable  flow  environment  within  the 
combustor.  There  is  a  need  to  understand  the 
response  of  the  shear  layer  to  externally  imposed 
acoustic  disturbances,  and  device  passive  and 
active  control  techniques  for  controlling  the 
behavior  of  this  flow  field.  Both  theoretical  and 
experimental  studies  are  being  carried  out  to 
understand  the  behavior  of  the  free  shear  layer 

within  dump  combustors. A  number  of  investi¬ 
gators  have  also  studied  numerically  the  flow 
fields  within  the  dump  combustors,  using  explicit 
time  marching  techniques,  with  and  without  viscous 


The  experimental  studies  reported  in  Ref.  1 
show  that  the  behavior  of  the  free  shear  layer  is 
sensitive  to  periodic  acoustic  disturbances 
imposed  at  the  downstream  boundary.  In  some  cases, 
it  has  been  shown  possible  to  drastically  alter 
the  size  of  the  recirculation  zone  using  disturb¬ 
ances.  There  is  a  need  to  systematically  study  the 
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significant  role  that  the  external  boundary 
conditions  and  the  related  acoustics  play  on  the 
behavior  of  the  free  shear  layer.  As  a  first  step 
towards  such  a  study,  the  flow  field  within  a  2-D 
dump  combustor  is  studied  through  the  numerical 
solution  of  the  unsteady,  compressible  Navier- 
Stokes  equations.  The  acoustics  characteristics  of 
the  inlet-combustor  system  is  altered  by  changing 
either  the  length  of  the  combustor,  or  altering 
the  downstream  boundary  condition,  by  partially 
blocking  the  downstream  boundary.  The  length  of 
the  inlet,  and  the  flow  Reynolds  number  are  held 
fixed  in  order  to  avoid  significant  changes  of  the 
thickness  of  the  shear  layer,  and  its  natural 
instability  characteristics. 

The  computed  flow  fields  are  analyzed  using 
iso-vorticity  plots  at  selected  time  levels.  The 
flow  properties  at  selected  locations  within  the 
shear  layer  are  also  analyzed  using  Fourier 
transform  techniques  to  identify  the  dominant 
frequencies.  The  shear  layer  downstream  of  the 
step  is  analyzed  using  classical,  linear  instabi¬ 
lity  analyses  to  identify  the  natural  frequency  of 
the  shear  layer.  It  is  found  that  the  boundary 
conditions  play  a  crucial  role  in  the  behavior  of 
the  unsteady  flow  within  the  combustor. 

NUMERICAL  FORMULATION 

In  Fig.  1,  the  three  configurations  being 
analyzed  are  shown.  The  typical  dimensions  of  the 
configurations  are  also  indicated,  normalized  with 
respect  to  the  step  height.  The  unsteady  flow 
within  the  combustor  is  governed  by  2-0  compress¬ 
ible  Navier-Stokes  equations,  and  is  likely  to  be 
turbulent.  Because  existing  algebraic  and  two 
equation  models  are  not  suitable  for  unsteady 
flows,  and  because  these  models  can  smear  out 
features  such  as  shear  layer  instability  and 
vortex  roll  up,  in  the  present  work  no  explicit 
turbulence  model  was  used.  An  algebraically 
generated,  stretched  Cartesian  coordinate  system 
was  used.  The  governing  equations,  which  are 
parabolic  with  respect  to  time,  were  integrated 
using  an  implicit  time  marching  procedure,  origi- 

O 

nally  devised  by  Beam  and  Warming.  This  procedure 
is  second  order  accurate  in  space,  and  may  be 
designed  to  be  either  first  or  second  order 
accurate  with  respect  to  time.  In  the  present 
work,  the  first  order  time  accuracy  option  was 
used.  This  procedure  has  been  previously  applied 
to  unsteady  external  flow  problems  with  good 
g 

success.  In  the  Appendix,  the  mathematical 
formulation  is  briefly  described. 

Since  numerical  solutions  are  sensitive  to 
grid  spacing,  a  variety  of  grid  sizes  were  exper¬ 
imented  with,  ranging  from  a  coarse  61  x  41  grid 
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system,  to  a  fine  151  x  121  grid  system.  As 
expected,  increased  grid  refinement  leads  to 
improved  resolution  of  features  such  as  the  shear 
layer  roll  up.  The  basic  characteristics  of  the 
flow  such  as  the  length  of  the  recirculation  zone 
at  a  given  time,  and  the  thickness  of  the  shear 
layer  are,  however,  insensitive  to  grid  spacing. 
Based  on  these  exploratory  studies,  in  the  calcu¬ 
lations  to  be  reported  here  it  was  decided  to  use 
the  151  X  121  grid  system.  Because  excessive  grid 
stretching  can  reduce  the  formal  spatial  accuracy 
of  the  solution  from  second  order  to  first  order, 
and  lead  to  a  loss  of  resolution  of  the  flow 
features,  a  uniform  grid  spacing  was  used  in  the 
studies  reported  here.  For  the  short  combustor 
configuration  shown,  this  is  equivalent  to  grid 
spacings  of  0.15  and  0.025  along  the  x-  and  z- 
directions  respectively. 

Numerical  viscosity  is  introduced  into  the 
solution  procedure  through  a  set  of  artificial 
viscosity  terms  as  explained  in  the  appendix.  In 
order  to  assess  the  effects  of  artificial  viscos¬ 
ity  on  the  solution,  calculations  were  performed 
for  values  of  the  artificial  viscosity  coefficient 
between  1  and  5.  The  flow  features  were  insensi¬ 
tive  to  this  coefficient,  within  this  range. 
Subsequently,  a  value  of  unity  was  used  for  the 
artifical  viscosity  coefficient,  in  the  calcula¬ 
tions  to  be  reported. 

BOUNDARY  CONDITIONS 

Because  the  governing  equations  are  parabolic 
with  respect  to  time,  the  proper  boundary  condi¬ 
tion  for  this  problem  is  the  specification  of 
velocity,  and  temperature  at  all  time  levels  at 
all  the  boundaries.  Unfortunately,  such  a  complete 
specification  of  the  boundary  conditions  is  seldom 
feasible,  and  rarely  available  from  experiments. 
Therefore,  the  following  set  of  approximate 
boundary  conditions  have  been  used. 

At  all  the  solid  walls,  the  no  slip  boundary 
condition  was  imposed.  Furthermore,  the  normal 
derivative  of  pressure  was  set  to  zero.  The 
temperature  at  the  solid  walls  was  evaluated  using 
adiabatic  assumptions.  The  density  at  the  wall  was 
subsequently  extracted  from  the  equation  of  state. 

At  the  inlet,  the  flow  was  assumed  to  be 
parallel  to  the  x-  axis,  and  the  velocity  profile 
was  assumed  to  be  a  "plug"  profile  (uniform 
everywhere  except  at  the  walls).  The  Mach  number 
at  the  inlet  was  chosen  to  be  0.2,  and  the  inflow 
density  was  assumed  to  be  unity.  Furthermore,  the 
derivative  of  pressure  along  the  x-  axis  was 
assumed  to  be  zero  at  the  inlet. 

At  the  outflow  boundary,  two  different 
treatments  are  needed  depending  on  whether  the 
downstream  boundary  is  partially  blocked,  or  is 
completely  open.  Configurations  1  and  2  shown  in 
Fig.  1  have  open  downstream  boundaries.  On  the 
open  boundary  the  pressure  was  prescribed,  while 
the  other  three  flow  variables  (density,  u  and  v 
components  of  velocity)  were  extrapolated  from  the 
interior.  The  portion  of  the  downstream  boundary 
that  is  blocked  was  treated  like  any  other  solid 
wall . 

The  above  boundary  conditions  are  considered 
"reflecting"  boundary  conditions.  That  is,  they 
allow  part  of  the  signals  attempting  to  leave  the 


computational  domain  to  be  reflected  back  at  the 
inlet  and  the  outflow  boundary.  In  a  real  world 
situation,  say  in  a  wind  tunnel,  no  such  reflec¬ 
tions  occur  although  there  may  be  other  sources  of 
disturbances  such  as  background  noise.  In  other 
cases,  say  within  a  ramjet  combustor,  the  inlet 
and  exit  may  be  choked.  Thus  the  above  boundary 
conditions  are  a  poor  substitute  for  what  takes 
place  near  the  boundaries  of  a  wind  tunnel  or  a 
ramjet  combustor.  As  mentioned  earlier,  removal  of 
false  reflections  from  the  boundaries  would 
require  prescription  of  velocity  and  temperature 
at  all  boundaries  at  all  time  levels. 

The  reflecting  boundary  conditions  described 
above  serve  one  very  useful  purpose,  however.  They 
provide  a  continuous,  small  supply  of  acoustic 
energy  at  a  frequency  which  is  characteristic  of 
the  dimensions  of  the  configuration.  These  waves 
can  interact  with  the  shear  layer,  and  lead  to 
amplification  of  shear  layer  instability,  shear 
layer  roll  up  and  other  interesting  phenomena.  For 
the  situation  where  the  pressure  is  fixed  at  the 
downstream  boundary,  and  has  zero  gradient  at  the 
upstream  boundary,  the  acoustic  mode  imposed  by 
the  reflecting  boundaries  corresponds  to  the 
quarter  wave  acoustic  mode.  Since  the  frequency 
of  this  mode  may  be  changed  by  changing  the  length 
of  the  configuration,  it  is  possible  to  study  the 
response  of  the  shear  layer  to  a  particular 
frequency,  simply  by  choosing  a  suitable  length  of 
the  configuration.  For  the  short  combustor  configu¬ 
ration  shown  in  Fig.  1,  the  quarter  wave  frequency 
happens  to  be  160  Hertz. 

RESULTS  AND  DISCUSSION 

Calculations  have  been  carried  out  for  the 
three  configurations  shown  in  Fig.  1.  In  all  the 
cases,  the  flow  was  started  impulsively,  assuming 
the  air  to  be^ stationary  within  the  configuration 
at  time  t  =  0  ,  except  at  the  inflow  boundary.  At 

subsequent  time  levels,  the  flow  velocity  within 
the  combustor  steadily  increases  as  the  stationary 
mass  of  air  is  replaced  by  the  high  speed  air- 
stream.  At  later  time  levels,  a  periodic  pattern 
is  established,  in  which  the  shear  layer  oscil¬ 
lates  in  resonance  with  a  frequency  that  is 
characteristic  of  the  combustor  length  for  config¬ 
urations  1  and  2,  and  the  cavity  length  for 
configuration  3.  During  a  given  cycle  a  part  of 
the  shear  layer  rolls  up  into  a  vortex,  and  is 
shed.  Depending  on  the  case  studied,  this  vortex 
may  merge  with  previously  shed  vortices  downstream 
at  a  subsequent  time. 

To  verify  the  observation  that  the  present 
boundary  conditions  lead  to  the  presence  of 
quarter  wave  acoustic  modes  for  configurations  1 
and  2,  a  large  number  of  calculations  were  done  on 
a  somewhat  coarser  grid,  for  a  number  of  config¬ 
urations.  The  length  of  the  combustor  was  paramet¬ 
rically  changed.  In  Fig.  2,  the  frequency  of 
flapping  of  the  shear  layer  is  plotted  as  a 
function  of  the  length  of  the  configuration.  The 
quarter  wave  acoustic  frequency  associated  with 
the  configuration,  given  by  4a/L  where  a  is  the 
speed  of  sound,  is  also  plotted.  It  is  seen  that 
these  two  quantities  match  well  for  the  entire 
range  of  combustor  lengths.  Similar  studies  have 
been  done  with  configuration  3,  which  indicate 
that  the  sh  >.sr  layer  oscillates  at  the  natural 
frequency  of  the  cavity. 
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Configuration  1:  For  a  detailed  discussion  of  the 
shear  layer  dynamics,  we  concentrate  on  configura¬ 
tion  1.  In  Fig.  3,  the  vorticity  contours  within 
the  dump  combustor  after  a  large  period  of  time 
are  shown,  at  regular  time  levels.  It  is  seen  that 
the  vorticity  contours  at  time  t  =  210  (normalized 
with  respect  to  step  height  and  inlet  speed  of 
sound)  are  identical  to  those  at  t  =  300.  The 
other  time  levels  (t  =  220,  t  =  310),  (t  =  230,  t 
=  320)  etc.  also  match.  That  is, the  shear  layer 
pattern  repeats  itself  once  every  90  non-dimen¬ 
sional  units  of  time.  This  exactly  equals  the 
quarter  wave  acoustic  frequency  of  configuration 
1,  which  has  a  length  equal  to  22.5. 

In  Fig.  3,  it  is  seen  that  the  shear  layer 
undergoes  considerable  lateral  oscillations  called 
flapping.  It  is  also  seen  that  the  shear  layer 
periodically  rolls  up  and  sheds  a  vortex,  for 
example  between  time  levels  220  and  290.  During 
the  time  interval  250  <  t  <  290,  this  vortex  pairs 
into  two  smaller  vortices  which  are  convected  out 
of  the  flow  domain. 

In  Fig.  4,  the  pressure  variation  at  a  point 
one  step  height  downstream,  and  in  the  middle  of 
the  shear  layer  is  plotted  as  a  function  of  time. 
It  is  seen  that  the  pressure  values  oscillate  at  a 
distinct  frequency,  but  the  amplitude  varies  from 
cycle  to  cycle.  The  calculations  were  carried  out 
for  approximately  12  cycles  of  oscillation  (approx¬ 
imately  1200  units  of  time,  24000  time  steps)  to 
ensure  that  the  phenomena  being  discussed  repeat 
themselves. 

In  order  to  understand  why  the  amplitude  of 
pressure  oscillations  vary  from  cycle  to  cycle, 
the  pressure  distribution  shown  within  4  was 
analyzed  using  Fast  Fourier  transform  techniques. 
In  Fig.  5,  the  Fourier  transform  of  the  pressure 
signals  at  six  different  locations  is  plotted. 
Note  that  the  distances  indicated  are  measured 
from  the  corner  of  the  step.  It  is  seen  that  the 
Fourier  transform  shows  pressure  peaks  at  four 
distinct  frequencies,  equal  to  160,  304,  416  and 
592  Hz.  The  fourth  frequency  appears  to  be  combi¬ 
nations  of  the  first  and  third  frequency  (160  + 
416  «  592).  The  160  Hz  peak  occurs  as  a  result  of 
the  quarter  wave  acoustic  mode. 

To  assess  the  reason  for  the  existence  of  the 
304  Hz  peak,  the  shear  layer  velocity  distribution 
just  downstream  of  the  boundary  was  analyzed  using 
classical  inviscid  shear  layer  stability  analysis. 
The  amplitude  of  the  streamwise  growth  of  disturb¬ 
ances  within  the  shear  layer  was  studied  as  a 
function  of  user  input  sinusoidal  temporal  varia¬ 
tion.  In  Fig.  6,  the  real  part  of  the  growth  rate 
is  plotted  as  a  function  of  the  input  frequency, 
for  a  velocity  profile  1.5  steps  downstream  of  the 
step.  It  is  seen  that  at  a  frequency  of  300,  the 
most  rapid  spatial  growth  of  the  linear  instabili¬ 
ty  waves  occur.  We  have  repeated  these  stability 
calculations  at  a  number  of  stations  in  the 
immediate  vicinity  of  the  step,  using  velocity 
profiles  selected  at  random  time  levels.  In  all 
cases,  the  analysis  indicated  that  at  300  Hz  the 
shear  layer  is  most  prone  to  instabilities. 

From  the  above  discussions  one  may  conclude 
that  the  shear  layer  within  the  configuration  1 
locks  onto  the  quarter  wave  acoustic  frequency  of 
the  shear  layer,  and  oscillates.  Shear  layer  roll 
up  and  shedding  of  vortices  also  occur  at  the 


quarter  wave  acoustic  frequency.  Secondary  oscil¬ 
lations  also  occur  within  the  shear  layer  at 
frequency  close  to  the  natural  frequency  of  roll 
up  of  the  shear  layer.  Some  of  the  higher  frequen¬ 
cies  observed  may  be  seen  to  be  combinations  of 
the  above  two  fundamental  frequencies  of  the  flow. 

Configuration  2:  This  configuration  differs  from 
configuration  1  only  in  the  length  of  the  com¬ 
bustor.  The  increased  length  of  the  combustor 
leads  to  a  somewhat  lower  quarter  wave  acoustic 
frequency,  equal  to  48  Hz.  In  Fig.  7,  the  vorti¬ 
city  contours  are  plotted  at  selected  time  inter¬ 
vals  for  the  period  200  <  t  <  300  .  The  natural 
frequency  of  roll  up  of  the  shear  layer,  computed 
using  velocity  profiles  just  downstream  of  the 
step  was  found  to  be  270  hz.  In  this  case  the 
shear  layer  still  locks  onto  the  quarter  wave 
acoustic  frequency  of  the  system  and  tends  to 
oscillate  very  slowly. 

In  Fig.  8,  the  Fourier  transform  of  the 
pressure  at  two  locations  within  the  shear  layer 
is  shown.  It  is  seen  that  two  distinct  peaks,  one 
at  the  system  acoustic  frequency  of  90  Hz,  and  the 
other  near  the  shear  layer  roll  up  frequency  of 
270  Hz  .  Because  the  time  step  was  chosen  to  be 
large  to  reduce  the  computer  time  requirements, 
higher  harmonics  or  combinations  of  these  fre¬ 
quencies  could  not  be  resolved  by  the  calcula¬ 
tions. 

Configuration  3:  This  configuration  differs  from 
the  previous  two  in  that  a  cavity  forms  between 
the  step  and  the  partially  blocked  downstream 
boundary.  The  cavity  has  a  distinct  acoustic 
frequency  of  its  own,  which  differs  from  the 
quarter  wave  acoustic  frequency  of  the  system 
based  on  the  length  of  the  configuration. 

Calculations  for  this  configuration  were 
done,  and  results  analyzed  using  techniques 
identical  to  those  employed  for  configurations  1 
and  2.  In  Fig.  9,  the  u-  component  of  velocity  and 
pressure  at  selected  locations  within  the  shear 
layer  is  plotted  as  a  function  of  time.  A  visual 
examination  of  this  near  sinusoidal  variation 
indicates  that  it  occurs  at  a  frequency  equal  to 
the  natural  acoustic  frequency  of  the  cavity.  In 
this  case,  the  velocity  profile  downstream  of  the 
step  tended  to  vary  rapidly  from  one  time  level  to 
the  other,  so  that  a  single  shear  layer  roll  up 
frequency  could  not  be  found.  A  Fourier  analysis 
of  the  velocity  and  pressure  variations  showed 
only  one  significant  peak,  corresponding  to  the 
natural  acoustic  frequency  of  the  cavity. 

CONCLUSIONS 

The  behavior  of  free  shear  layers  within  2-D 
dump  combustors  has  been  analyzed  using  numerical 
solution  of  time-dependent  Navier-Stokes  equa¬ 
tions.  The  computed  flow  fields  have  been  analyzed 
using  vorticity  contour  plots.  Fast  Fourier 
transform  of  the  pressure  fluctuations,  and 
linear,  inviscid  stability  analysis  of  the  shear 
layer.  The  following  conclusions  may  be  drawn, 
based  on  the  present  study. 

1)  Boundary  conditions  play  a  crucial  role  in  the 
behavior  of  the  flow  within  the  combustor.  For 
example,  when  the  dev/nstream  boundary  is  blocked 
the  shear  layer  tends  to  oscillate  at  the  cavity 
frequency,  rather  than  the  system  frequency. 
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Changing  the  length  of  the  configuration  leads  to 
large  variations  In  the  system  response,  when 
monitored  as  a  function  of  time. 

2)  The  shear  layer  tends  to  lock  onto  the  natural 
acoustic  frequency  of  the  system  (or  of  the 
cavity).  The  flapping  motion  of  the  shear  layer 
was  accompanied  by  vortex  shedding,  and  pairing  In 
the  cases  studied. 

3)  Fourier  spectra  of  the  pressure  fluctuations 
within  the  system  show  a  second,  somewhat  weaker 
peak,  at  the  shear  layer  Instability  frequency. 
Some  of  the  higher  frequencies  found  In  the 
Fourier  spectra  appear  to  be  combinations  or 
multiples  of  these  two  basic  frequencies. 
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according  to  the  following  one-to-one  relation¬ 
ship: 


C=C{x)  ;  n=n(y)  :  T=t  (1) 

The  Jacobian  of  transformation  J  Is  given  by 


J  = 


^x'^y 


(2) 


and  the  metrics  of  transformation  are  given  by 


(3) 


Standard  central  differences  were  used  to 
compute  quantities  such  as  x  ,  x  etc.  which  In 
return  were  used  to  compute  quantities  such  as 
,  t  etc.  At  the  solid  surface  and  the  inflow- 
/oull^Flow  boundaries,  three-point  one-sided  dif¬ 
ferences  were  used  to  compute  the  metrics. 


In  the  (C,n,i)  coordinate  system,  the  two- 
dimensional,  unsteady  Navier-Stokes  equations  may 
be  written  as 
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where 

q  =  J  ^{p,pu,pv,e}  (5) 

The  quantities  F,  G,  R  and  S  are  given  by 


F  =  F  )  /  J 
6  =  (hy  G)  /  0 
»  *  R)  /  J 
S  *  (hy  S)  /  J 


F= 


PU2 
pu  +p 


puv 

Lu(e+p) J 


G= 


pv 

puu 

pv  +p 

v(e+p) 


(6) 


(7) 


The  quantities  R^  and  S^  represent  the 
dissipation  of  energy  ^due  to  work  done  by  the 
viscous  stresses,  and  heat  conduction  along  the  x- 
and  y-  directions  respectively.  The  viscous 
stresses  t  ,  t  and  t  were  related  to  the 
velocity  giYdlents  through^ Stokes'  hypothesis.  As 
mentioned  earlier  the  objective  of  the  present 
work  Is  to  determine  the  onset  of  laminar  shear 
layer  Instability,  and  no  explicit  turbulence 
model  has  been  used  In  the  computations  to  be 
presented  here. 


9.  Sankar,  L.N.  and  Tang,  W.,  "Numerical  Solu¬ 
tion  of  Unsteady  Viscous  Flow  past  Rotor 
Sections,"  AIAA  paper  85-0129. 

APPENDIX 

All  the  calculations  were  performed  In  a 
stretched  Cartesian  coordinate  system  (C.n,T) 
which  Is  linked  to  the  physical  coordinate  system 


DISCRETIZATION  AND 
APPROXIMATE  FACTORIZATION 

Since  the  governing  equations  are  coupled  to 
each  other  and  are  highly  nonlinear,  a  stable, 
efficient  solution  procedure  Is  required  for 
solving  them.  In  the  present  work,  the  Beam- 
Warming  algorithm  was  used  with  some  modifica¬ 
tions.  The  viscous  terms  were  explicitly  evaluated 
using  Information  available  at  earlier  steps. 
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since  the  mathematical  and  numerical  formulation 
of  the  Beam-Warming  algorithm  are  well  known,  only 
a  brief  description  of  the  solution  scheme  is 
given  here. 

The  governing  equations  are  written  at  a 
computational  node  (i,j)  in  the  following  finite 
difference  form: 

+  6-F  .  ,  G  =  S-R  +  6  S  -  Ec  D 

e  +  5n  t  n  E 

where  for  example,  the  term  a  is  the  stand¬ 
ard  two  point  central  difference  formula  given  by 
(F.  ,  -  Fs.j)"  /2.  The  quantity  D  is  the  artifi¬ 
cial^  dissipation  term  discussed  in  the  next 
section. 

The  highly  nonlinear  terms  F  and  G  at  the 
time  level  (n-*-!)  were  expanded  by  a  Taylor  series 
about  a  previous  time  level  n  as  shown  below: 

=  f"  +  [DF/Dq]"  -  q") 

g"^^  =  g"  +  [OG/Oq]"  (q""^  -  q")  (9) 

Here  the  quantities  OF/Dq  and  OG/Oq  are  4x4 
matrices  which  are  the  Jacobians  of  the  flux 
vectors  F  and  G  with  respect  to  q. 

In  order  to  allow  large  values  of  the  expli¬ 
cit  dissipation  coefficient  Cc  to  be  used  with  out 
instability,  and  to  allow  thf  viscous  terms  to  be 
treated  explicitly,  the  following  implicit  dissi¬ 
pation  terms  were  added  to  the  left  side  of  the 
difference  equation  (8): 

J  (q"*^  -  q”)  (10) 

The  coefficient  e.  was  taken  to  be  three 
times  the  explicit  dissipation  coefficient  c,  .  A 
range  of  Cp  values  between  3  and  5  were  usSd  in 
the  calculations  reported  here. 

Equation  (8)  may  be  written  after  the  addi¬ 
tion  of  the  artificial  implicit  dissipation  terms 
given  by  Equation  (10),  in  the  following  operator 
form: 

[I+4t6^  {DF/Dq}+Jt«^ {OG/Oq} 

-eiJ'^{«55+«„,,)J](«l"''^-q")=  R"  (11) 

where, 

R"  =  -dt(4,F  +  «  G)"  +4t(6,.R 
5  n  5 

4^$)"  -dtEg  O"  (12) 

The  left-hand  side  operator  of  Equation  (11) 
was  approximately  factored  into  two  smaller 

operators,  leading  to  the  following  final  form: 

[I+4t4^{0F/0q}-ejdtj‘^4^^J] 

[I+dt4  (OG/Oq) 

n 

-Cjitj‘^4^^J](q'’''^  -  q")=  r"  (13) 

Equation  (13)  may  be  solved  through  the 

inversion  of  two  block  tridiagonal  matrix 
equations,  one  corresponding  to  the  (-  direction 
and  the  other  corresponding  the  n~  direction.  In 


order  to  keep  the  flow  solver  simple,  the  boundary 
conditions  on  all  the  boundaries  were  explicitly 
updated  after  the  interior  points  had  been  updated 
using  Equation  (13). 

TREATMENT  OF  THE  EXPLICIT 
DISSIPATION  TERMS 

In  the  earlier  Euler  and  Navier-Stokes 
equations,  researchers  used  the  following  form  of 
the  artificial  dissipation  term: 

0"  -Jt  (J,)"  (14) 

This  term  is  formally  of  the  order  of  the  fourth 
power  of  the  grid  spacing  in  the  physical  plane, 
and  is  not  expected  to  reduce  the  overall  accuracy 
of  the  solution  technique.  This  form  was  found  to 
give  nonphysical  overshoots  in  the  vicinity  of 
rapid  flow  gradients  such  as  shocks.  It  was  found 
that  second  order  artificial  dissipation  terms  did 
not  exhibit  a  similar  overshoot,  but  led  to  highly 
inaccurate  solutions. 

A  solution  to  the  problem  of  overshoots  was 
proposed  by  Jameson.  In  his  approach,  the  dissi¬ 
pation  term  was  written  as  a  combination  of  second 
and  fourth  order  dissipation  terms.  A  sensor, 
based  on  the  second  derivative  of  pressure  turned 
on  the  second  order  dissipation  in  the  vicinity  of 
shocks,  and  suppressed  the  fourth  order  dissi¬ 
pation  term.  Away  from  the  rapid  gradients,  the 
fourth  order  dissipation  form  was  used.  Jameson's 
approach  was  implemented  in  the  following  study  as 
f o1 1 ows . 

The  term  D  was  written  as 

0=4tJ’^[4  -  AC14„+ 

*■  n{A(l-Cl)4nnn  n 

4^{B(1-C2)4^^^-B  C24^}](Jq)"  (15) 

The  coefficients  C2and  C,  are  proportional  to 
the  second  derivatives  of pressure,  and  are 
defined  such  that  it  will  of  significant  value  (of 
the  order  of  unity)  only  near  rapid  gradients  such 
as  shock  waves.  Elsewhere,  these  coefficients  are 
of  the  order  zero,  and  the  expression  given  in 
equation  (15)  leads  to  a  fourth  order  error  in  the 
solution.  The  coefficients  A  and  B  are  propor¬ 
tional  to  the  wave  speed  in  the  x-  and  y-  direc¬ 
tions  and  provide  an  upwind  flavor  to  the  present 
scheme . 
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Figure  3.  Sample  Vortlclty  Contours  at  Selected  Time  Levels  for  Configuration  1. 
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Figure  4.  Pressure  Time  History  at  a  Point  Within 
the  Shear  Layer  (x/h  =  1.05,  y/h  =  0.0)  for 
Configuration  1. 
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Figure  5.  Fourier  Spectrum  of  the  Pressure 
Fluctuations  within  the  Shear  Layer  for 
Configuration  1. 


frequency,  Hz 

•Figure  6.  Amplitude  of  Streamwise  Disturbance 
Growth  within  the  Shear  Layer  as  a  Function  of 
Frequency  for  Configuration  1. 
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Figure  7.  Sample  Vorticity  Contours  at  Selected 
Time  Levels  for  Configuration  2. 


Figure  8.  Fourier  Transform  of  Pressure 
Fluctuations  within  the  Shear  Layer  for 
Configuration  2. 


Figure  9.  Velocity  and  Pressure  Fluctuations  at 
Selected  Points  within  the  Shear  Layer  for 
Configuration  3. 
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ABSTRACT 

The  Issue  of  enhancing  mixing  between 
parallel,  supersonic  streams  is  numerically 
investigated.  An  explicit  time  marching  scheme  that 
is  second  order  accurate  in  time  and  fourth  order 
accurate  in  space  is  used  to  study  this  probleta 
Small  amplitude  velocity  disturbances  at  selected 
frequencies  are  imposed  over  an  otherwise  steady 
flow  at  the  juncture  of  the  two  streams  to  prrxTX^e 
mixing.  It  is  found  that  disturbances  are  selectively 
amplified  at  certain  frequencies,  while  disturbances 
at  other  frequencies  are  rapidly  damped  out  In 
studies  where  the  relative  Mach  number  the 
disturbances  relative  to  one  of  the  streams  is  high, 
shocklets  were  found  to  form  on  one  or  both  sides 
of  the  shear  layers.  In  such  a  situation,  the  relative 
Mach  numbers  of  the  eddies  were  different  in 
coordinate  systems  attached  to  the  upper  and  the 
lower  streams. 

INTRODUCTION 

Aircraft  engine  and  missile  manufacturers 
are  presently  interested  In  a  class  of  propulsion 
systems  called  SCRAMJET  engines.  In  these 
systems  the  supersonic  airstream  captured  at  the 
inlet  is  slowed  down  to  modest  supersonic  speeds 
through  a  series  of  shock  waves  prior  to  entering 
the  combustion  chamber.  Here  the  airstream  is 
allowed  to  mix  and  react  with  a  parallel  stream  of 
fuel  or  partially  burnt  fuel/air  mixture.  For  efficient 
performance  of  these  systems,  it  is  necessary  that 
the  fuel  arxl  air  streams  mix  with  each  other  as 
rapidly  as  possible,  over  a  foiriy  short  distance. 

Unfortunately,  supersonic  free  shear  layers 
which  form  at  the  juncture  of  the  air  and  fuel 
streams  terxJs  to  grow  very  slowly  (Ref.  1-3] 
compared  to  their  subsonic  counterparts.  Alternate 
mechanisms  such  as  normal  injection  of  fuel  into 
the  airstream  will  likely  increase  mixing,  but  at  the 
expense  of  significant  total  pressure  losses. 
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Therefore,  there  is  some  interest  in  the  use  of  active 
and  passive  control  techniques  which  will  promote 
mixing. 

PREVIOUS  WORK 

A  comprehensive  discussion  of  recent 
experimental,  numerical  and  analytical  studies  on 
the  behavior  of  subsonic  and  supersonic  shear 
layers  has  been  done  by  Dimotakis  [Ref.  4],  Here, 
only  a  small  subset  of  existing  work,  doseiy  related 
to  the  present  numerical  studies,  is  reviewed. 

Experimental  Studies:  Chinzei  et  al.  [Ref.1]  have 
experimentally  studied  the  growth  rate  of  planar 
shear  layers,  using  Schiieren  techniques,  and  total 
pressure  probes.  They  found  organized  vortical 
structures  to  exist  in  such  flows,  in  a  manner  similar 
to  subsonic  shear  layers.  Perhaps  the  best  known 
experimental  work  on  supersonic  planar  shear 
layers  is  that  done  by  Papamoschou  and  Roshko 
(ref.  2,3],  for  a  variety  of  gases  and  flow  conditions 
on  either  side  of  the  shear  layer.  They  showed  that 
the  convective  Mach  numbw  of  the  eddies  is  a 
significant  parameter  governing  the  growth  rate  of 
supersonic  shear  layers.  Papamoschou  also 
performed  stability  analyses  of  infinitely  thin  shear 
layers  (vortex  sheets)  to  link  the  growth  rate  of  the 
shear  layer  (compart  to  that  of  an  incompressible 
shear  layer)  to  the  convective  Mach  number,  and 
derived  do^  form  expressions  for  the  convective 
Mach  number  as  a  fuiiction  of  flow  conditions  on 
either  side.  The  idea  of  convective  Mach  number 
itself  b,  of  course,  not  new,  and  has  been 
previously  derived  by  Bogdanoff  [Ref.  5].  In  a  later 
study  [Ref.  6],  Papamoschou  found  that  the 
measured  corrvective  Mach  number  of  the  eddies 
matches  the  analytical  predictions  only  when  the 
convective  Mach  number  is  low  and  subsonic.  He 
attributed  this  discrepancy  to  the  foct  that  the 
traditional  derivations  for  the  convective  Mach 
number  assume  the  total  pressure  on  either  side  of 
the  shear  layer  to  be  equal.  In  cases  where  the 
convective  Mach  number  is  high,  shocklets  can 
form  arxf  lead  to  different  amounts  of  total  pressure 
losses  on  etther  side  of  the  shear  layer. 
Papamoschou  also  studied  modifications  to  the 
traling  edge  of  the  splitter  plate  which  initiaily 
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In  order  to  understand  the  behavior  of 
supersonic  free  shear  layers  at  low  convective 
Mach  numbers,  we  investigate  its  response  to 
arbitrary,  user-specified  acoustic  disturbances  over 
a  broad  range  of  frequencies.  Sinusoidally  varying 
velocity  disturbances  at  a  number  of  frequencies 
are  introduced  at  the  initial,  laminar  mixing  region  of 
the  shear  layer.  These  disturbances  grow  with  time 
as  they  are  convected  downstream  and  eventually 
lead  to  well  organized  vortical  structures.  The 
objective  of  this  work  is  then  to  study  how  the 
disturbances  over  the  entire  spectrum  of 
frequencies  behave  as  they  are  convected 
downstream,  and  to  speculate  on  mechanisms  by 
which  energy  is  transferred  from  high  frequencies 
to  low  frequencies  and  vice  versa. 

To  study  behavior  of  the  shear  layer  at  very 
high  convective  Mach  numbers,  we  use  vorticity 
and  pressure  contour  plots  at  a  number  of  time 
levels  to  track  the  velocity  of  the  dominant  eddies 
arxj  compute  the  relative  Mach  number  of  these 
eddies  in  a  coordinate  system  attached  to  either  the 
faster  stream  or  the  slower  stream.  If  supercritical 
Mach  numbers  arise  relative  to  either  stream,  then 
the  resultant  pressure  field  is  examined  for  the 
occurrence  of  shock  waves,  expansion  waves  and 
their  effects  on  the  shear  layer  growth. 

The  2-0  compressible  Navier-Stokes 
equations  in  a  strong  conservation  form  are 
numerically  solved,  using  a  modified  MacCormack 
scheme  that  is  second  order  accurate  in  time,  and 
fourth  order  accurate  in  space.  This  scheme  Is 
suitable  for  studying  phenomena  such  as 
propagation  of  acoustic  waves,  boundary  layer 
instabOity,  and  shear  layer  instability  and  has  been 
prevlousiy  used  by  several  authors  [Ref.  22-24].  The 
flow  field  is  assumed  to  be  laminar. 

NUMERICAL  FORMULATION 

The  2-0,  laminar,  unsteady,  compressible 
flow  Is  governed  by  the  Navier-Stokes  equations 
which  may  be  formally  written  as: 

qj  +  Fj^  +  Gy-R^  +  Sy 

where  F  and  G  are  inviscid  flux  terms,  while  R  and  S 
are  the  viscous  stress  terms. 


In  this  work  the  above  equation  was  solved 
using  a  splitting  approach.  That  is,  the  solution  was 
advanced  from  one  time  level  'n'  to  the  next  Cn+2') 
through  the  following  sequence  of  operations: 

^  ^  *-y  *T(v  ^  ^(v  *-y  *-x 

where,  for  example,  the  Lx  operator 
involves  solution  of  the  following  1  -0  equation: 

qt  +  F^  -  0 

This  1-0  equation  was  solved  through  the 
following  predictor-corrector  sequence, 
recommended  by  Bayliss  et  al  [Ref. 22]: 


Predictor  Step: 

q,*.  q,"  iAt/(a^x)  [7  Fj  -  8  Fj.,  +  F^gl" 


Corrector  Step: 

q,"+^.(q,%q,")/2+At/(12ax)[7Fi.8Fj^^  +Fi^2l* 

In  the  above  equations,  the  j-  index  has 
been  suppressed  for  darfty. 

When  the  above  equation  is  applied  at 
nodes  dose  to  the  left  and  right  side  boundary,  a 
fourth  order  accurate  extrapolation  procedure  wa| 
used  to  extrapolate  the  flux  vectors  F  and  F 
needed  at  nodes  outside  the  computational 
domain. 

The  Ly  operator  requires  solution  of  the 
equation 

qj  +  Gy  «  0 

using  a  simflar  approach. 

The  operators  and  correspond  to  numerical 
solution  of  1-0  equations  such  as 

qt-R^-O 
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equivalent  to  applying  a  slip  boundary  condition  at 
these  walls.  In  the  case  of  free,  unconfined  shear 
layers,  non-reflective  boundary  conditions  are 
needed  at  these  lateral  boundaries.  In  the  present 
work,  setting  the  derivative  of  the  normal 
component  of  velocity  to  zero,  rather  than  the 
velo^  Itself  to  zero  was  found  to  minimize 
reflections. 

The  flow  properties  at  the  inflow  were 
specified  everywhere  in  the  flow  field  as  the  initial 
conditions  for  the  problem.  The  Navier-Stokes 
solver  was  then  advanced  for  several  non- 
dimensional  units  of  time,  until  a  fully  developed 
shear  layer  with  a  modest  streamwise  growth  was 
established. 

Once  a  steady  state  shear  flow  was 
achieved,  forced  excitation  of  the  shear  layer 
began.  This  was  achieved  by  prescribing  the 
normal  (v-)  component  of  velocity  over  the  entire 
inflow  boundary  to  behave  as  follows: 

v(y.x-0.t)  -  Z  An  f(y)  sin^^nt+tf  „) 

Here  the  summation  shown  is  over  all  the 
excitation  frequencies;  is  the  amplitude  of 
disturbance,  Un  is  the  frequency  of  disturbance  and 
0  P  is  the  associated  phase  an^e.  The  function  f(y) 
determines  the  variation  of  the  perturbation  velocity 
across  the  shear  layer.  Both  a  Gaussian  distribution 
and  a  constant  magnitude  distribution  were 
attempted.  The  results  to  be  presented  here 
correspoTKl  to  f(y)  equal  to  unity. 

In  the  present  work  6  frequencies  were 
used,  with  zero  phase  difference  between  the 
individual  components.  The  quantity  Ap  was  2%  of 
the  reference  velocity  U^.  The  non-dimensional 
frequencies  w-  were  10,  20,  30,  40.  50  and  60 
respectively.  Obviously,  a  linear  stabflity  analysis 
could  have  been  used  to  pick  the  frequencies  that 
are  related  to  the  most  unstable  frequency.  But  the 
intent  here  was  to  impose  excitation  at  several 
frequencies  on  the  shear  layer,  at  the  Inflow 
boundary  and  determine  which  frequencies  are 
selectiveiy  amplified,  and  to  determine  what 
happens  to  the  energy  content  at  the  higher 
frequencies  at  subsequent  time  levels. 

Figure  2  shows  the  vorticity  contours  at  a 
randomly  selected  time  level.  It  Is  seen  that  the 


vorticity  field  at  the  immediate  downstream 
bound^  is  rich  in  structure  showing  large 
gradients  in  the  streamwise  as  well  as  normal 
directions.  At  large  distances  downstream, 
however,  only  a  single  row  of  eddies  at  well  defined 
distances  are  seen. 

Because  the  formation  and  motion  of 
vortices  (or  eddies)  give  rise  to  a  rapidly  varying 
pressure  field  which  moves  with  the  eddies  some 
useful  information  about  the  energy  content  at  the 
shear  layer  distributed  over  the  various  frequencies 
may  be  obtained  by  computing  the  Fourier 
transform  of  the  pressure  field  at  a  number  of  points 
within  the  shear  layer.  In  figure  3,  the  Fourier 
spectrum  of  the  pressure  field  is  plotted  at  6  x- 
locations  within  the  computational  field,  at  y/S  «  0. 
The  following  trend  is  seen.  Near  the  inflow 
bourxlary,  the  Fourier  spectrum  ^^hows  a  near 
uniform  distribution  over  the  entire  frequency  range. 
At  downstream  locations,  the  higher  frequency 
content  begins  to  gradually  decrease.  The  low 
frequency  components  at  norvdimensionai 
frequencies  10  and  20  show  a  rapid  increase 
Initiaily,  but  reaches  asymptotically  constant  values. 
Rgure  4  contains  the  same  information  as  figure  3, 
except  it  shows  the  changes  in  energy  content  as  a 
function  of  downstream  distance.  It  is  seen  that 
Fourier  coefficient  associated  with  non-dimensional 
frequency  60  reduces  to  30%  of  initial  value  100  5 
downstream,  whereas  the  low  frequency 
component  triples  in  magnitude  and  reaches  its 
limit  value  125  units  downstream.  An  examination 
of  the  vorticity  contour  plot  (figure  2)  shows  a 
number  of  sniall  eddies  at  the  inflow  boundary, 
which  rapidly  merge  into  a  single,  large  vortex.  This 
merging  appears  to  be  the  mechanism  responsible 
for  the  decrease  in  the  energy  content  at  high 
frequencies,  arxf  the  corresponding  increase  at  the 
lower  frequencies.  It  is  interesting  to  note  that  the 
two  lowest  frequency  components  (corresponding 
to  non-dimensional  frequencies  10  and  20)  maintain 
their  energy  levels  once  they  reach  their  limit 
values,  with  no  further  transfer  of  energy  from  the 
Up  -  20  waves  to  thew^*  10  waves.  This  may  be 
due  to  the  fact  that  the  phase  difference  in  the 
forcing  function  corresponding  to  these  two  waves 
was  zero. 
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In  figures  7  and  8.  the  pressure  and  velocity 
contours  are  plotted  at  a  randomly  chosen  time 
level.  In  this  case,  from  an  inspection  of  the  vorticity 
profiles  at  adjacent  time  levels,  the  eddies  appear  to 
travel  at  sutetantially  lower  speed  than  the  upper 
stream.  During  the  early  stages  of  eddy  formation 
and  motion,  shocidets  occur  both  on  the  upper  and 
lower  sides  of  the  shear  layer.  As  the  eddies 
accelerate  and  reach  low  subsonic  Mach  numbers 
relative  to  the  upper  stream,  the  shocidets  on  the 
upper  side  of  the  shear  layer  disappear.  The 
shMidets  on  the  lower  side  continue  to  travel  with 
the  eddies,  with  no  reduction  in  their  strength. 

Case  3:  Case  1  was  repeated,  by  artificiaiiy 
reducing  the  shear  layer  vorticity  thickness  by  a 
factor  of  IS,  keeping  ^1  (Mher  dimensions  such  as 
the  grid  size,  domain  length  arxl  width  constant  In 
figures  9  and  10,  the  pressure  and  vorticity 
contours  are  plotted.  Again,  in  figure  10,  the  solid 
and  dotted  contours  corresporxJ  to  low  and  high 
pressure  levels  respectively.  Shocidets  are  evident 
on  either  side  of  the  shear  layer,  although  they  are 
weak  because  the  reduced  shear  layer  thickness 
leads  to  small  and  thin  eddies,  compared  to  case  1. 

Case  4:  As  a  final  exercise.  Case  1  was  repeated, 
with  forced  excitation  of  the  normal  velocity  at  the 
inflow  boundary  over  multiple  frequencies,  ranging 
from  10  to  60.  The  amplitude  of  the  irxlividual 
components  was  0.02  times  the  upperstream 
veio^.  In  figure  11  the  Fourier  spec^m  of  the 
pressure  field  at  several  x-  locations  are  plotted.  A 
gradual  migration  of  energy  levels  from  the  higher 
frequency  to  the  lower  frequencies  is  evident  as  in 
the  case  of  the  subsonic  convective  Mach  number 
case.  Figure  12  shows  how  the  two  high  frequency 
components  decay  following  a  brief  initial  growth  as 
eddies  are  convected  downstream.  The  low 
frequency  components  at  frequencies  10 , 20  and 
30  initially  grow  rapidly,  but  reach  asymptotic 
values. 

The  response  of  the  shear  layer  to  multiple 
frequencies  is  strikingly  simflar  to  that  in  the  earlier 
study  for  subsonic  convective  Mach  numbers, 
shotm  In  figures  3  and  4. 

CONCLUDING  REMARKS 

The  stability  arxl  growth  characteristics  of 
supersonic  free  shear  layers  were  studied  by 


exciting  the  shear  layer  at  the  upstream  boundary 
with  small  amplitude  normal  velocity  disturbances. 
The  foliowing  observations  were  made: 

a)  In  the  case  of  shear  layers  at  subsonic 
and  supersonic  convective  Mach  numbers,  the 
Imposition  of  acoustic  disturbances  over  a  large 
range  of  frequencies  lead  to  the  trartsfer  of  this 
energy  from  the  high  frequencies  to  the  low 
frequencies,  as  the  flow  progressed  from  the 
upstream  bourxlary  to  the  downstream  boundary. 
The  energy  content  at  the  lowermost  frequencies 
rapidly  reached  asymptotic  values  following  which 
eddies  in  the  shear  layer  were  convected 
downstream  with  no  further  alteration  in  their 
structure. 

b)  In  the  case  of  shear  layers  at  a 
supersonic  convective  Mach  number,  situations 
were  found  where  the  convective  Mach  number 
relative  to  the  faster  stream  is  low.  This  leads  to  a 
situation  where  shocidets  arose  only  on  the  lower 
side  of  the  shear  layer.  Corxjitions  were  also  found 
where  the  convective  Mach  number  relative  to  both 
the  streams  is  high,  leading  to  shockiets  on  either 
side.  These  calculations  demonstrate  the  same 
features  experimentally  observed  by  Papamoschou 
[Ref.  3]  and  discussed  based  on  total  pressure 
arguments  by  DinrK>takis  (Ref  4]. 
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Figure  1 .  Computational  domain 


Figure  4.  Variation  of  Fourier  spectrum  as  a  function  of 
streamwise  location;  McsO.2.  Mi»4.0,  M2>2.3. 
ai/a2-2.3,  5-1/15. 


Figure  2.  Vorticity  Contours  for  a  shear  layer  excited  at 
multiple  frequencies;  McaO.2,  Mi-4.0,  M2-2.3, 
ai/a?— 2.3, 5-1/15. 
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Figure  3.  Fourier  spectrum  of  pressure  field  at  selected 
x-locations  within  the  shear  layer;  Mc-0.2. 
Ml— 4.0,  M2— 2.3,  ai/a2— 2.3, 5—1/15. 
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Rgureii.  Fourier  spectrum  of  pressure  field  at 
selected  x-locations  within  the  shear  layer; 
Mca1.2,  Mla6.0,  M2s3.6,  ai/ajal.O,  &b1/15. 


Figure  12.  Variation  of  Fourier  spectrum  as  a  function  of 
streamwise  location;  Mc>l.2,  Mis6.0, 
M2a3.6,  ai/a2>1.0,  5b1/1S. 
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Numerical  Simulation  of  the  Growth  of  Instabilities  in 
Supersonic  Free  Shear  Layers 


W.  Tang,*  N.  M.  Komerath.f  and  L.  N.  Sankart 
Georgia  Institute  of  Technology,  Atlanta,  Georgia 


The  behavior  of  the  initial  region  of  a  snpcnonic  plane  shear  layer  is  analyzed  ibrongh  namerical  solution  of 
the  two-dimensional  Navier-Stokes  equations,  as  weii  as  the  three-dimensional  eqnations  under  the  finite-span 
assumption.  A  modified  MacCormadi  scheme  that  is  fourth-order  accurate  in  space  and  second-order  in  time 
is  employed.  Small  amplitnde  oscillations  in  the  normal  vdodiy  arc  fonnd  to  grow  as  they  coavcct  downstream, 
and  eventually  lead  to  organized  vortical  structures.  Normal  velocity  disturbanecs  are  found  to  be  more  efficient 
than  streamwise  or  spanwisc  disturbances.  The  growth  rate  of  these  disturbances,  as  well  as  the  intensity  of 
velocity  fluctuations,  are  fonnd  to  decrease  at  the  convective  Mach  number  of  the  shear  layer  increases.  The 
Mach  number  of  the  vortical  structures  with  respect  to  the  faster  stream  is  found  to  be  considerably  less  than 
the  theoretical  value  of  the  convective  Mach  anmher. 


Nomenclature 

a\ 

=  upper  stream  speed  of  sound 

3  lower  stream  speed  of  sound 

F,C 

=  flux  vectors 

Afvwta 

-convective  Mach  number  of  present  result 

K 

-convective  Mach  number 

g 

-  vector  of  conserved  variables 

R,S 

-  diffusion  vectors 

Uc 

-convection  speed  of  vortices 

C/| 

-  upper  stream  inflow  velodty 

t/j 

=  lower  stream  inflow  velocity 

X 

-  coordinate  in  streamwise  direction 

y 

-  coordinate  in  normal  direction 

Introdaction 

IR-BREATHING  engines  designed  for  high-flight  Mach 
numbers  require  supersonic  combustion  for  efficient  op¬ 
eration.  The  shock  losses  associated  with  deceleration  to  low 
Mach  numbers  require  that  the  mixing  of  fuel  and  air,  and  the 
heat  release,  must  occur,  in  supersonic  flows.  For  the  same 
reason,  it  is  desirable  to  mix  the  fuel  and  air  using  coflowing 
streams.  In  such  conngurations,  the  mixing  must  occur  across 
the  shear  layer  formed  between  the  streams.  The  length  ard 
weight  of  the  engine,  and  the  efficiency  of  heat  release,  depend 
on  the  rapidity  of  this  mixing  process.  Most  current  concepts 
for  supersonic-combustion  ramjets  thus  employ  mixing-lim¬ 
ited  heat  release.  The  mixing  across  a  shear  layer  between  two 
streams  depends  on  the  rate  of  mass  and  momentum  transfer 
across  the  layer,  and,  hence,  can  be  described  using  the 
“growth”  or  “spreading”  rate  of  the  shear  layer.  Unfortu¬ 
nately,  shear  layers  separating  supersonic  streams  are  known 
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to  grow  much  more  slowly  than  corresponding  subsonic  shear 
layers. 

One  long-term  objective  of  supersonic  shear  layer  research, 
therefore,  is  to  devise  methods  of  increasing  the  mixing  be¬ 
tween  supersonic  streams  by  enhancing  the  shear  layer  growth 
rate.  Recent  success  in  greatly  modifying  subsonic  shear  layers 
has  resulted  in  the  advancement  of  a  variety  of  schemes  for 
achieving  similar  increases  in  supersonic  shear  layers.  The 
variety  of  such  possibilities  far  exceeds  the  resources  available 
for  experiment^  exploration  of  each.  Instead,  a  better  ap¬ 
proach  appears  to  be  to  develop  reliable  numerical  models  and 
solution  methods  that  can  then  be  used  to  perform  the  explo¬ 
ration,  and  to  identify  promising  approaches  and  the  appro¬ 
priate  values  of  parameters  required.  This  is  the  motivation 
behind  the  research  described  in  this  paper. 

Prevhw  Work 

Chinzei  et  al.'  conducted  experiments  on  planar  shear  layer 
configurations  and  studied  the  growth  rate.  Papamoschnu^ 
conduaed  similar  experiments,  using  a  variety  of  gases  and 
flow  conditions,  and  showed  that  the  results  could  be  scaled 
using  the  convective  Mach  number  of  the  dominant  eddies  in 
the  shear  layer.  These  results  showed  that  the  growth  rate  of 
supersonic  shear  layers  is  typically  less  than  one-third  the 
growth  rate  of  incompressible  shear  layers  for  convective 
Mach  numbers  greater  than  unity. 

Passive  and  active  control  techniques  have  been  studied  by 
otl  r  researchers.  These  techniques  are  generally  based  on  the 
principle  that  if  vortidty  is  introduced  into  the  shear  layer,  it 
will  increase  the  level  of  fluctuation  and,  therefore,  promote 
mixing  and  growth.  Guirguis^  and  Drummond  and  Mukunda* 
studied  the  effect  of  a  bluff  body  placed  in  the  middle  of  the 
shear  layer.  Kumar  et  al.’  considered  the  effects  of  vortidty 
produced  by  a  pulsating  shock  wave  on  the  growth  character¬ 
istics  of  the  shear  layer.  Ragab  and  Wu*  have  developed 
calculations  based  on  stability  theory  to  predict  the  response 
of  supersonic  shear  layers.  Recently,’  they  have  also  developed 
computations  of  the  response  of  planar  wakes  and  shear  layers 
similar  to  those  in  experimental  splitter  plate  configurations. 

Scope  of  Prcecol  Paper 

In  the  work  presented  here,  the  behavior  of  a  planar  free 
shear  layer  is  studied,  using  two  numerical  techniques  for 
solution  of  the  Navier-Stokes  equations.  The  effects  of  active 
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control  strategies  are  investigated.  Sinusoidal  variations  in  the 
velocity  are  introduced  at  the  upstream  boundary.  The  subse* 
quent  response  of  the  shear  layer  to  these  disturbances  is 
studied.  Streamwise,  normal,  and  spanwise  disturbances  are 
considered  as  suitable  candidates  for  promoting  mixing. 

At  present,  the  problem  is  assumed  to  be  nominally  two-di¬ 
mensional.  Some  calculations  have  been  performed  with 
three-dimensional  layers  under  the  infinite  sweep  assumption. 
It  is  recognized  that  the  later  development  of  the  shear  layer 
may  be  strongly  influenced  by  three-dimensional  effects. 
However,  there  is  no  reason  to  believe  that  the  initial  region 
should  be  anything  other  than  two-dimensional.  The  available 
experimental  flow  visualizations,  performed  with  spanwise-in- 
tegrating  techniques  such  as  schlieren  and  shadowgraphy, 
clearly  show  structures  that  would  have  been  totally  smear^ 
out  if  the  flowflelds  had  been  significantly  three-dimensional. 

The  present  calculations  are  for  laminar  :hear  layers,  and 
no  turbulence  model  is  used.  Turbulence  models  inherently 
bring  additional  uncertainty  into  the  physical  interpretation  of 
the  observed  behavior  of  the  flowfleld,  though  they  are  cer¬ 
tainly  necessary  to  obtain  quantitative  accuracy.  The  lack  of 
such  a  model  restricts  the  applicability  of  these  results  to  the 
initial  region  of  the  shear  layer. 

The  initial  velocity  profile  used  is  a  step  change  in  velocity 
at  the  slip  line  between  the  two  streams.  Thus,  the  results 
obtained  will  not  correspond  to  experimental  results  from 
splitter-plate  configurations,  since  there  is  no  boundary  layer 
and  no  embedded  region  of  initially  subsonic  flow. 

Within  the  above  limitations,  the  present  work  aims  to 
study  the  behavior  of  the  initial  region  of  a  shear  layer,  and  to 
explore  the  effects  of  various  forms  of  excitation. 


Problem  Statement 

The  shear  layer  configuration  is  shown  in  Fig.  1.  Two 
uniform,  parallel  supersonic  streams  of  different  Mach  num¬ 
bers  are  released  at  the  left-hand  boundary.  All  properties  are 
known  at  this  boundary.  The  upper  and  lower  boundaries  of 
the  computational  domain  are  assumed  to  be  hard  walls  across 
which  no  disturbances  can  escape.  There  is  no  boundary  layer 
at  these  walls,  and  slip  conditions  are  used.  At  the  downstream 
boundary,  the  flow  and  all  disturbances  are  allowed  to  escape, 
and  no  disturbances  are  allowed  to  propagate  back. 

To  study  shear  layer  behavior,  the  static  pressures  are  equal¬ 
ized  across  the  splitter  plate,  so  that  there  are  no  strong  shocks 
in  the  flow.  Some  waves  and  their  reflections  from  the  wall  do 
occur,  but  these  are  quite  weak. 

The  flow  is  assumed  to  be  nonreacting,  and  the  ratio  of 
specific  heats  was  assumed  to  be  constant  for  both  streams. 
The  species  above  and  below  the  shear  layer  were  assumed  to 
have  the  same  molecular  weight. 

Maibcinetical  Fonnuielion:  Fourth  Order  MecCormack  Scbcaae 

g,+Fj+Gy^Rx+S,  (1) 


Py  =  Uy=Py=Ty=0,  V=0 


Py  =  “  y  =  Py  =  Ty  =  0,  V  =  0 


Fig.  1  Boundary  coaditioaa  for  tupcnonk  free  shear  layer. 


Here,  F  and  G  are  the  inviscid  flux  terms  and  account  for 
the  transport  of  mass,  moment,  and  energy  and  for  the  influ¬ 
ence  of  pressure.  The  terms  R  and  S  are  the  viscous  stress 
terms.  The  above  equation  is  parabolic  with  respect  to  time, 
and  may  be  solved  using  a  variety  of  stable  time  marching 
schemes.  For  two-dimensional  flows,  there  are  four  equations. 
In  the  case  of  three-dimensional  flows  subject  to  infinite- 
sweep  assumption,  there  are  five  equations,  the  additional 
equation  corresponding  to  the  conservation  of  spanwise 
momentum. 

In  this  work,  the  above  equation  was  solved  using  a  splitting 
approach;  that  is,  the  solution  was  advanced  from  one  time 
level  ‘n '  to  the  next  ‘n +2' ,  through  the  following  sequence  of 
operations: 

q''*^  =  {LyLyLxyL,J.,J.rJ.yLj,)q'’  (2) 

where,  for  example,  the  L,  operator  involves  solution  of  the 
following  equation: 

q,+F,=0  (3) 

This  one-dimensional  equation  was  solved  through  the  fol¬ 
lowing  predictor-corrector  sequence,  recommended  by  Bayliss 
et  al.*: 

Predictor  Step: 

9o-  =  P-j-— ^  ^F,j-iF,.ij+F,.2j  (4) 


Corrector  Step: 

‘  ^  (ffu" + •)  +  7^  - 1 J  +  F, .  ij  ' 

(5) 

When  the  above  equations  are  applied  at  nodes  close  to  the 
left-  and  right-side  boundary,  a  fourth-order  accurate  extrapo¬ 
lation  procedure  was  used  to  extrapolate  the  flux  vectors  Fand 
F'  needed  at  nodes  outside  the  computational  domain. 

The  Ly  operator  requires  solution  of  the  equation 

q,  +  C,=0  (6) 

using  a  similar  approach. 

The  operators  L„  and  Ly,  correspond  to  numerical  solution 
of  one-dimensional  equations  such  as 

q,-R.=0  (7) 


The  above  equation  was  solved  through  the  following  two-step 
sequence: 


The  viscous  terms  are  thus  updated  only  to  second-order 
accuracy  in  space.  It  may  be  shown  that  the  above  scheme  has 
very  little  artificial  dissipation  inherent  in  it,  and  is  fourth-or¬ 
der  accurate  in  space,  as  far  as  the  inviscid  part  is  concerned. 
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Boundary  Conditions 

As  suted  above,  all  flow  properties  are  prescribed  at  the 
upstream  boundary  for  both  streams,  including  any  imposed 
perturbations.  At  the  downstream  boundary,  the  flow  is  as¬ 
sumed  to  remain  fully  supersonic  for  the  small-amplitude 
perturbations  encountered  in  this  work,  so  that  the  properties 
may  be  extrapolated  from  the  interior.  Alternatively,  the  gov¬ 
erning  equations  themselves  may  be  applied  if  the  streamwise 
diffusion  terms  R,  are  suppress^  at  the  downstream  nodes. 

At  the  lateral  boundaries,  the  flow  is  assumed  to  be  con¬ 
fined  by  smooth,  parallel  walls.  Slip  boundary  conditions  were 
used  to  avoid  the  compression  effects  that  would  be  caused  by 
boundary  layers.  The  walls  were  considered  adiabatic,  and  the 
normal  derivatives  of  density  and  pressure  were  set  to  zero. 

Results  and  Discussion 

Normalization 

Velocities  were  normalized  using  the  speed  of  sound  in  the 
upper  stream,  which  thus  became  unity.  The  Reynolds  number 
based  on  the  speed  of  sound  was  chosen  to  be  1000.  A  221  x  61 
uniform  grid  was  used,  with  spacing  of  0.01  in  each  direction. 
Thus,  the  length  of  the  domain  L  was  2.2.  The  time  step  was 
taken  as  0.001  (0.0005  for  each  half-step).  The  calculations 
were  started  with  step  velocity  profiles  at  the  upstream 
boundary,  and  allowed  to  proceed  until  a  steady  state  was 
reached  asymptotically.  This  usually  took  6(X)  time  steps.  The 
results  at  this  stage  were  stored,  and  the  code  was  restarted 
with  an  imposed  sinusoidal  velocity  disturbance  of  amplitude 
2*70  of  the  velocity  of  the  upper  stream.  The  calculations  were 
then  run  until  several  cycles  of  the  disturbance  had  been 
completed,  and  the  initial  effects  had  been  convened  away 
through  the  downstream  boundary. 

Convective  Mach  Number 

The  cases  run  have  been  summarized  in  Table  1 .  Because  the 
supersonic  shear  flow  problem  involves  several  parameters  (at 
least  five  on  either  side  of  the  shear  layer),  nondimensional 
groupings  are  sought  to  express  observed  effects.  Following 
the  practice  of  Papamoschou,^  the  convective  Mach  number 
was  used  here.  For  the  problem  studied  here,  the  convective 
Mach  number  reduces  to 

A/c  =  (f/.-C/c)/ai 

where 

Ue  =  (Oi  C/j  +  a2f/i)/(0|  +  Oj) 

The  values  of  Me  calculated  by  this  formula  are  tabulated.  A 
physical  interpretation  of  the  convective  Mach  number  is  that 
it  is  the  Mach  number  of  the  dominant  large-scale  vortical 
structures  with  respect  to  either  stream.  According  to  the 
formula  given  above,  this  Mach  number  is  the  same  with 
respect  to  either  of  the  streams.  An  attempt  was  made,  as 
discussed  later,  to  determine  the  convection  speed  of  the  vorti¬ 
cal  structures  seen  in  the  computational  flowfield,  and  to 
determine  relative  Mach  numbers  from  them.  The  Mach  num¬ 
bers  so  determined,  with  respect  to  the  upper,  high-speed 
stream,  are  also  tabulated.  It  is  seen  that  there  is  a  consider¬ 
able  discrepancy.  This  is  not  surprising,  and,  in  fact,  even  in 
subsequent  experiments  by  Papamoschou,’  similar  effects  ap¬ 
pear  to  have  been  observed. 


Table  1  Cases  presented 


Case 

Ml 

M2 

Oi 

£/j 

Me 

A/vocte* 

1 

4.0 

2.3 

4.00 

3.51 

0.20 

0.2 

2 

4.0 

2.0 

4.00 

3.05 

0.38 

0.2 

3 

4.0 

1.3 

4.00 

1.98 

0.80 

0.2 

4 

5.0 

1.3 

S.OO 

1.98 

1.20 

0.6 

Fomution  of  Vortical  Stmcturcs 

Figure  2a  shows  the  contours  of  vorticity  in  the  shear  layer, 
calculated  for  case  1,  with  no  disturbance  superposed.  It  is 
seen  that  the  shear  layer  grows  quickly  at  the  very  beginning, 
and  then  takes  on  a  smooth  profile  which  grows  very  little 
thereafter.  It  should  be  remembered  that  in  this  calculation 
there  is  no  imposed  turbulence  model.  Figure  2b  shows  the 
effect  of  imposing  a  sinusoidal  2*7s  normal  velocity  distur¬ 
bance  at  the  inflow  boundary.  Distinct  centers  of  vorticity  are 
seen  to  develop  and  be  convected  downstream.  The  shear  layer 
edge  now  penetrates  considerably  further  into  both  streams. 
Careful  examination  of  the  contours  shows  considerable 
asymmetry  and  distortion  as  the  structures  proceed  down¬ 
stream.  The  computational  domain  in  this  calculation  does 
not  extend  far  enough  for  these  disturbances  to  grow  into  the 
nonlinear  regime,  and  hence  no  “roll-up”  can  be  expected 
here.  The  effects  of  six  cycles  of  the  imposed  disturbance  can 
be  seen,  with  the  sixth  just  leaving  the  computational  domain. 

Figures  3a  and  3b  show  the  corresponding  vonicity  con¬ 
tours  for  case  2,  where  the  theoretical  value  of  convective 
Mach  number  is  nearly  twice  that  of  case  1 .  The  growth  rate 


x/H 

Fig.  2  Vortkity  conloars  for  case  1,  M|  •  4.0,  M:  •  2.3,  Mt »  0.2.; 
a)  wiiboot  dtotarbaacc,  b)  wiib  dblurbaacc  in  ibe  aormsd  direction. 
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F1f.3  Vortidty  eontoan  forcaae  2,  Mr»0.38: 

a)  wllboni  dMarbance,  b)  with  dbiarbance  In  ibe  normal  dircctioa. 
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x/H 

Fig.  4  Vortldty  contours  for  c«*«  3,  A/i«4.0,  1.3,  Afe»0.8; 

a)  without  disturbance,  b)  with  disturbance  in  the  normal  direction. 


.3 


-.1 


-.2- 

-.3  •  I  I  r  I . .  '  ■  I 
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Fig.  5  Vortkity  contours  for  case  4,  Afi*5.0,  4fi«U,  Mt^l.2; 
a)  without  disturbance,  b)  with  dbtnrbance  in  the  normal  direction. 


.3- 

.2- 

.1 


-.1  - 
-.2- 

'^  0  2  4  .6  .8  1.0  1.2  1.4  1.6  1.8  2.0  2.2 

x/H 

Fig.  6  Effect  of  kreamwbe  disturbances  for  case  1. 
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Fig.  7  Effect  of  spanwise  disturbances  for  case  1. 
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Fig.  8  Compulation  of  convective  speed  based  on  the  temporal  evo- 
lation  of  vortkity  contours  for  case  2. 


appears  to  be  less,  as  expected  from  the  experimental  observa¬ 
tions  of  the  effect  of  A#,.  This  effect  is  seen  further  in  Figs.  4 
and  5,  where  the  convective  Mach  number,  according  to  the 
formula,  is  0.8  and  1.2,  respectively. 

Distnrbancc  Type 

Other  kinds  of  disturbance  were  tried;  specifically,  distur¬ 
bances  in  velocity  along  the  streamwise  direction  for  the  two- 
dimensional  shear  layer  and  disturbances  along  the  spanwise 
direction  using  a  three-dimensional  model  with  the  infinite 
sweep  assumption.  The  results  are  shown  in  Figs.  6  and  7, 
respectively,  for  the  case  where  the  convective  Mach  number  is 
predicted  to  be  0.2.  In  each  case,  the  disturbance  amplitude  is 
the  same.  It  is  seen  that  these  types  of  disturbances  are  less 
efficient  than  the  normal  velocity  disturbance. 


Convective  Speed  Based  on  Evolution  of  Vortkity  Contours 
One  common  way  of  calculating  the  convective  Mach  num¬ 
ber  is  to  track  the  downstrebm  convection  of  the  vorticity 
contours  as  a  function  of  time.  The  results  of  this  effort  are 
shown  in  Fig.  8  for  case  2.  Similar  computations  were  per¬ 
formed  for  all  of  the  cases,  and  the  results  are  shown  in  Table 
I.  It  is  seen  that  as  the  theoretical  value  of  the  convective 
Mach  number  rises,  the  measured  convective  Mach  number  of 
the  structures  with  respect  to  the  upper  stream  vary  much  less. 
Of  course,  this  means  that  the  convective  Mach  number  is  no 
longer  the  same  when  compared  to  the  two  streams.  The 
structures  appear  to  move  at  a  speed  that  is  close  to  that  of  the 
upper,  high-speed  stream.  This  is  similar  to  the  more  recent 
observations  of  Papamoschou,*  where  the  structures  were 
found,  for  many  cases,  to  move  at  speeds  close  to  that  of  one 
or  the  other  stream. 
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Growth  Rate  Eahancemciit 

The  thickness  of  the  shear  layer  was  computed  from  the 
velocity  profiles  across  the  shear  layer.  These  profiles  are 
shown  later  in  the  paper,  where  they  are  used  to  examine  the 
numerical  accuracy  of  the  results.  Figure  9  shows  the  shear 
layer  thickness  for  unperturbed  and  perturbed  cases  for  case  I. 
The  disturbed  case  shows  a  significantly  greater  rate  of 
growth,  except  near  the  downstream  boundary.  However,  the 
increase  in  growth  rate  is  only  on  the  order  of  10-15*70. 

The  frequency  of  the  imposed  fluctuations  was  chosen  such 
that  about  six  vortical  structures  could  be  seen  in  the  computa¬ 
tional  domain  at  one  time.  Thus,  the  actual  frequency  used 
was  higher  than  the  frequency  of  maximum  ampliHcation 
predicted  by  linear  stability  analysis.  Use  of  the  preferred 
frequency  would  have  required  a  much  larger  computational 
domain  in  the  x  direction  to  capture  an  adequate  number  of 
vortices,  and  would  have  increased  the  computational  effort 
required  by  an  order  of  magnitude. 


0.2  0.4  0.6  0.8  1.0  1.2  1.4  1.6  1.8  2.0  2.2 
x/H 

Fig.  10  Effect  of  coovcetive  Macb  naabcr  oo  flnctaatioo  la  the 
shear  layer;  a)  M,  -0.2;  b)  M,  >0.30;  aad  c)  Mt  -O.M. 


The  intensity  of  fluctuation  in  the  shear  layer  is  plotted  as  a 
function  of  downstream  distance  for  three  cases  with  different 
theoretical  convective  Mach  numbers  in  Figs.  lOa-lOc.  The 
quantity  measured  was  the  root-mean-squared  fluctuation  of 
the  {/-component  of  velocity  about  the  mean,  normalized  by 
the  mean  velocity.  It  is  to  be  noted  that  this  is  not  the  turbu¬ 
lence  intensity,  since  no  attempt  has  been  made  to  model  the 
turbulence.  The  fluctuations  have  their  origin  in  the  imposed 
disturbance,  though  they  may  have  been  selectively  ampliHed 
by  energy  exchange  with  the  shear  layer.  It  is  seen  that  the 
intensity  of  fluctuations  decreases  rapidly  with  increasing  con¬ 
vective  Mach  number.  It  also  appears  that  the  intensity 
quickly  reaches  an  asymptotic  value  and  does  not  increase 
further.  In  fact,  for  the  higher  values  of  convective  Mach 
number,  the  intensity  appears  to  peak  and  then  decrease  grad¬ 
ually  thereafter.  This  decay  in  the  intensity  with  high  values  of 
Me  has  been  predicted  by  other  researchers  using  linear  stabil¬ 
ity  analysis. 

Studies  of  DIscrelization  Error 

All  the  above  results  were  generated  using  the  fourth-order 
MacCormack  scheme.  The  influence  of  the  accuracy  of  the 
computation  scheme  was  studied  by  comparing  results  ob¬ 
tain^  using  a  second-order  MacCormack  scheme  to  those 
obtained  with  the  fourth-order  scheme.  Velocity  profiles 
across  the  shear  layer  were  used  to  examine  the  results.  Figure 
11  shows  the  comparison  for  case  1,  however,  with  a  111x31 
grid.  The  profile  was  obtained  at  the  station  10<7o  of  the 
domain  length  downstream  of  the  origin.  The  results  are  seen 
to  be  quite  similar.  The  agreement  is  close  for  proHles  at  25 
and  50*7o  downstream.  However,  at  x/L=0.75,  differences 
can  be  clearly  seen.  The  fourth-order  scheme  is  seen  to  resolve 
spatial  details  better,  as  expected.  The  difference  decreased 
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Fig.  II  Velocity  profiles  across  the  shear  layer:  comparison  of  sec¬ 
ond-  and  fonrtb-order  MacCormack  Kbeaies  for  case  I  with  III  x  31 
grid. 
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Fig.  12  Velocity  profiies  ecroes  the  shear  iayer:  comparisoB  of  sec¬ 
ond*  and  fourth-order  MacCormack  schemes  for  case  1  with  221  x  61 
grid. 


further  downstream.  It  is  concluded  from  these  that  the 
second-order  scheme  already  shows  reasonable  accuracy,  and 
that  the  fourth-order  results  are  quite  accurate. 

The  effea  of  grid  size  was  checked  by  comparing  the  above 
results  with  those  obtained  with  the  221  x  61  grid,  as  shown  in 
Fig.  12.  In  this  case,  the  difference  between  the  second-  and 
fourth-order  codes  is  negligible  for  all  stations.  Thus,  it  is 
concluded  that  with  the  221x61  grid,  the  results  are  not 
sensitive  to  grid  size. 

Compniationai  Resources 

All  calculations  reported  here  were  performed  on  the  Cray 
X-MP  at  the  Pittsburgh  Supercomputing  Center.  The  CPU 
time  per  time  step  per  grid  node  was  10  |ts  for  the  fourth-order 
MacCormack  scheme,  using  the  221  x  61  grid. 

Conclusions 

A  numerical  study  of  the  behavior  of  planar  supersonic 
shear  layers  has  been  performed.  Different  numerical  schemes 
have  been  tested.  Techniques  for  enhancing  the  growth  rate  of 
the  shear  layer  have  been  investigated.  The  following  are 
noted  from  the  results: 

l)The  fourth-order  MacCormack  scheme  accurately  simu¬ 
lates  the  evolution  of  the  imposed  disturbances.  The  results 
are  essentially  independent  of  the  spatial  order  of  the  scheme, 
with  a  221  X  61  grid  with  the  grid  parameters  used. 


2)  A  perturbation  of  2V*  of  the  mean  flow  velocity,  imposed 
at  the  upstream  boundary  in  the  normal  direaion,  produces  a 
larger  growth  of  the  shear  layer  than  an  equal  amount  of 
perturbation  imposed  in  the  streamwise  or  spanwise  direc¬ 
tions. 

3)  Imposed  sinusoidal  disturbances  in  the  normal  velocity 
upstream  lead  to  the  formation  and  growth  of  vortical  struc¬ 
tures.  The  shear  layer  thickness  grows  rapidly  at  first  and  then 
the  growth  rate  decreases  asymptotically. 

4)  The  root-mean-square  fluctuation  level  in  the  streamwise 
velocity  and  the  shear  layer  growth  rate  decrease  with  increas¬ 
ing  values  of  the  theoretical  convective  Mach  number  of  the 
shear  layer. 

5)  The  vortical  structures  are  found  to  move  at  different 
Mach  numbers  relative  to  the  upper  and  lower  stream,  and  the 
relative  Mach  number  appears  to  be  smaller  relative  to  the 
stream  with  the  higher  Mach  number. 
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Abstract 

The  temporal  stability  and  growth  characteristics  of  three- 
dimensional  supersonic  shear  layers  are  numerically  in¬ 
vestigated.  An  explicit  time-marching  scheme  that  is  second- 
order  accurate  in  time  and  fourth-order  accurate  in  space  is 
used  to  study  this  problem.  The  shear  layer  is  excited  by 
insubility  waves  computed  from  a  linear  stability  analysis  and 
random  initial  disturbances.  At  low  convective  Mach  num¬ 
bers,  organized  vortical  structures  develop  both  for  the  ran¬ 
dom  disturbance  and  the  modal  disturbance  cases.  At  super¬ 
sonic  convective  Mach  numbers,  vortical  structures  devdop 
initially  but  are  not  sustained  in  time.  Temporal  growth  of 
disturbances  is  found  to  be  a  strong  function  of  the  convective 
Mach  number. 


The  Ly,  and  operators  involve  solutions  of  the  one-di¬ 
mensional  equation  such  as 

(3) 


This  one-dimensional  equation  is  solved  through  the  following 
predictor-corrector  sequence,  recommended  by  Bayliss  and 
Maestrello.* 

Predictor  step: 

q;  -  qr  -  J'  W 


Corrector  step: 

q<  -  2 


Dr 

12  Dx 


17F,  -  8F,* 


(5) 


Coatents 

An  improved  understanding  of  factors  that  contribute  to 
supersonic  shear-layer  growth  is  necessary  for  design  of  active 
and  passive  control  techniques  to  enlmce  the  mixing  of 
airstreams  and  fuel  streams,  and  for  the  design  of  effldent, 
compact  SCRAMJET  engines.  It  has  been  observed'*^  that,  in 
supersonic  shear  layers,  organized  vortical  structures  exist  in  a 
manner  similar  to  subsonic  shear  layers.  However,  as  the 
convective  Mach  number  inaeases,  the  streamwise  shear-layer 
growth  rate  is  found  to  drop  to  about  309t  of  that  of  an 
incompressible  flow.* 

In  the  past,  Tang  et  al.’  used  a  fourth-order  MacCormack 
scheme  to  study  temporal  and  spatial  growth  of  two-dimen¬ 
sional  thin  shear  layers  at  very  early  stages  of  laminar  mixing, 
and  studied  the  effects  of  convective  Mach  number  as  well  as 
streamwise,  spanwise,  and  cross-stream  velocity  disturbances 
on  the  shear  layer  growth.  It  was  demonstrated  that  the 
growth  rate  of  the  shear  layer  decreased  with  increasing  con¬ 
vective  Mach  number.  In  this  work,  three-dimensional,  tem¬ 
porarily  growing  mixing  layers  have  been  studied.  The  study 
focuses  on  the  effects  of  instability  waves  computed  using  a 
linear  stability  analysis  and  random  initial  disturbances  on  a 
temporarily  evolving  shear  layer. 

liie  three-dimensional,  laminar,  unsteady,  compressible 
flow  is  governed  by  the  Navier-Stokes  equatioiu,  which  may 
be  formally  written  in  a  strong  conservation  form 

Qt  +  Fj,  +  Gy  +  Hi  •  Rx  +  Sy  +  T,  (I) 

where  F,  G,  and  H  are  invisdd  flux  terms,  and  R.  S,  and  Fare 
the  viscous  stress  terms.  Equation  (1)  was  solved  using  an 
operator  splitting  approach  and  a  MacCormack-type  finite 
difference  scheme: 

q"  *  *  —  ‘jf ’ll 'jnlvJ -iJ »>i,f -aXif- *  (2) 
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The  viscous  operators  L„Ly„  and  L„  are  integrated  simi¬ 
larly  with  the  exception  that  the  viscous  stress  terms  are  differ¬ 
ence  in  space  with  the  second-order-accurate  central  differ¬ 
ence  scheme.  The  overall  numerical  scheme  is  fourth-order 
accurate  in  space  and  second-order  accurate  in  time  as  far  as 
the  invisdd  part  is  concerned. 

Since  the  temporal  development  mixing  layers  are  studied, 
periodic  boundary  conditions  in  the  streamwise  and  spanwise 
directions  and  slip  boundary  conditions  in  the  cross-stream 
direcdon  are  applied. 

The  compuutional  domain  is  a  rectangular  channel  that 
extends  in  stream-  and  spanwise  directions  over  one  wave¬ 
length  of  the  longest  disturbance  wave  predicted  by  linear 
stability  analysis  for  a  given  convective  Mach  number.  In  the 
cross-stream  direction,  it  extends  from  -7.3  to  7.3  times  the 
vortidty  thickness.  The  computational  domain  is  discretized 
with  a  66  X  34  X  121  uniformly  spaced  grid  along  the  stream- 
wise,  spanwise,  and  cross-stream  directions,  respectively.  The 
Reynolds  number  is  based  on  the  vortidty  thickness  and  ' 


Fig.  I  Modal  coergy  growth  of  the  RMtt  ooMaMe  owdas  In  ahtar 
laytn  dWorhed  by  hMUMMly  waves. 
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ranges  between  3  x  10^  and  6  x  10>.  The  mean  velocity  across 
the  cross-stream  direction  is  given  by  a  hyperbolic  tangent 
profile.  The  convective  Mach  number  is  deHned  as  A#,  - 
(U\  -  Ua/iCx  •¥  cH,  where  {/  is  the  magnitude  of  the  mean 
velocity,  c  is  the  speed  of  sound,  and  subscriptt  1  and  2  refer 
to  the  upper  and  lower  streams,  respectively. 

The  modal  kinetic  energy  content  of  the  flowfleld  is  de- 
Hned  as 

-  jlioi*  +  w*  -(•  *rwT  di  (6) 

where  « ,  v ,  and  w  are  the  two-dimensional  Fourier  transforms 
of  the  velocity  field  on  the  plane  that  spans  in  the  streamwise 
and  spanwise  directions.  The  integration  is  in  the  CTOSS-stream 
direction.  The  superscript  *  denotes  the  complex  conjugate. 

lastaHmy  Waves  Sapwpeeed  on  Mean  Flow 

We  have  first  superposed  three-dimensional  waves  onto  the 
mean  flow  and  monitored  temporal  evolution  of  the  flowfield 
and  the  modal  energy  growth.  These  disturbance  waves  are  the 
most  unsuble  waves  predicted  by  the  linear  stability  analysis^ 
at  a  given  convective  Mach  number  and  are  given  for  the 


Fig.  2  Sganwies  vortMty  esnlanw  at  aMmon  la  a  raadaady  dla- 
tarbid  maar  hqwr,  M,  -  1.2. 


Fig.  3  Average  peitnrbadoa  kiactk  energy  growta  in  randoady  dis- 
larbed  shear  layers. 


velocity  components,  density,  and  temperature  in  the  follow¬ 
ing  form: 

d(x,y.z)  *  AD(z)  exp  /(ax  +  0y)  (7) 

where  D{z)  is  the  eigenfunction,  a  and  0  are  the  wave  num¬ 
bers,  and  ^  is  the  magnitude,  which  is  set  to  O.OlSA/r. 

We  have  studied  the  flowfields  for  convective  Mach  num¬ 
bers  of  0.2,  0.7,  and  1.2.  The  corresponding  wavelengths  of 
the  most  unstable  waves,  which  are  predicted  by  the  linear 
subility  analysis  code,  are  a  «  0.41  and  /3  -  0  for  A/r  >  0.2, 
a  >0.3  and/3  0.3  for  A/r  >  7,  and  a  >  0.14 and  d  «  0.07  for 
Mt  >  1.2.  The  temporal  growth  of  modal  kinetic  energy  asso¬ 
ciated  with  the  most  unstable  modes  is  shown  in  Fig.  1 . 

Raadow  Distarhaaces  Sagsrgestd  oa  Mcaa  Flow 

Next,  we  superposed  a  random  initial  disturbance  field  onto 
the  mean  flowfleld.  The  random  disturbance  field  was  gener¬ 
ated  using  a  random  number  generator,  and  its  magnitude  at 
any  point  in  the  flowfleld  was  restricted  to  be  less  than 
0.03Mr-  These  disturbances  were  confined  to  regions  of  signif¬ 
icant  vorticity  in  the  shear  layer,  where  lu(y)l  sO.lSMf. 
Computations  were  performed  for  conveaive  Mach  number  i; 
of  0.2,  0.4,  0.6,  and  1.2. 

In  all  flow  cases,  random  organized  vortical  structures  were 
observed  in  the  perturbation  velocity  field.  However,  at  a 
higher  convective  Mach  number  of  1.2,  organized  struaures 
tended  to  die  out  in  time  (Fig.  2).  The  temporal  growth  of 
average  perturbation  kinetic  energy  for  the  cases  studied  is 
given  in  Fig.  3. 

The  following  is  concluded: 

1)  The  temporal  growth  rate  in  perturbation  kinetic  energy 
decreases  with  increasing  convective  Mach  numbers  for  both 
modal  and  random  disturbances. 

2)  At  supersonic  convective  Mach  numbers  (A/r  «  1.2),  the 
growth  of  three-dimensional  structures  were  also  found  to  be 
unsustainable  in  both  the  random  and  modal  disturbance 
cases. 
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ABSIRACT 

A  tedmiffie  for  eonqmting  2>D  and  3-D  nmteady, 
ineoiqvessible  vuoous  flows  is  described.  Tbis  procedure  allows 
fourth  or  sixth  order  accuracy  in  space,  and  second  order  accuracy 
in  time,  and  is  general  enough  to  be  used  to  study  external  and 
intenal  flow  problems.  Algorithmic  details  of  the  procedure  are 
described.  A  saaqrle  application  of  this  algorithm  to  3-D  steady 
triscous  flow  over  a  boefy  of  revolution  is  given.  The  final  paper 
wfll  give  additional  nummeal  results  for  external  flow  past  a  body 
of  revolution  as  an  angle  of  attadt,  and  flow  over  2-D  airfoil 
geometries. 


INTRODUenON 

The  problem  of  accurately  conqiuting  unsteady  three 
dhnensiooal  incompressible  flows  is  one  of  great  interest  to 
rescardiers  working  in  the  fields  of  bioengineering  and  aeronace 
industiy.  The  flow  within  biftircated  oannels,  flow  in  blood 
vessels  with  con^liant  waUs,  flow  over  conmiex  submarine 
shapes,  flow  past  underwater  propellers,  and  flows  within 
tnrbomadiinery,  ducts  and  piunps  are  inherently  three- 
dimensional,  a^  unsteady.  These  conjSgurations  can  not  be 
poperiy  analyzed  unless  a  computaocual  c^iabiliw  that 
emblements  eryeiimental  methods  exists  to  analjRK  such  flows. 

Tedmiqnes  for  efBdent  and  accurate  prediction  of  3'D, 
fncon^ressibl^  unsteady  flows  are  also  necessary  fnr  the  first 
princ^  based  prediction  of  incompresrible  tubulem  flows 
tfaroi^  direct  numerical  simulations  or  large  eddy  simulaticns. 
The  ude  of  sndi  tools  is  one  of  the  princ^  reasons  that  first 
princ^  based  prediction  of  turbulent  flows  past  and  throng 
oompiaz  oonfifutations  have  not  been  extensively  attenqrted  to 


Finally,  flierc  has  been  a  significant  interest  within  the 
leisardi  eoinmunity  during  the  past  several  years  in  the  use 
active  and  nassive  contnd  tedmimes  to  improve  flow  quality  and 
reduce  undesirable  dienomena.  There  is  crasiderable  interest  in 
the  use  of  control  aeviett  to  promote  ntixing  in  air/fliel  streams, 
or  bot/cold  streams,  to  ensure  uniform  inflow  at  tlw  compressor 
fact,  and  to  prevem  flow  separation  within  diflbers,  ano  highly 
curved  ducts.  There  is  also  interest  in  the  use  df  control  devices  to 
reduce  noise,  vibration,  and  flutter.  These  control  devices  and 


control  laws  can  not  be  perfected  until  a  o^ability  fat  accurate 
and  efficient  prediction  of  unsteady  incmqiressible  flows  is 
available. 


A  SURVEY  OF  EXISITNO  SOLUTION  PROCEDURES 


Interest  In  the  drvekpment  of  numerical  solution 
procedures  for  the  prediction  of  two-  and  three-dimensional 
unsteady  Navier-Stokes  equations  arose  in  the  mid  sixties  and 
ear^  seventies  with  the  availability  of  high  qpeed  digital  com¬ 
puters  of  the  IBM  370  and  the  CDC  6600  dasL  During  the  past 
two  decades,  evolution  of  these  techniques  has  been  Aided  by  the 
availabiliqr  of  faster  and  dieaper  scalar/vector  computer 
architectures,  the  availabiliqr  of  effident  sohiuon  algorithms,  and 
the  develobirat  of  automated  grid  generation  tediniques 
capable  of  dividing  conqrlex  flow  regions  into  snioothly  vaiying 
cells.  At  the  presmit  time,  this  field  of  computational  fluid 
namia  it  a  weD  developml  sdence,  documented  in  text  books 
(Ref.  1-3],  conference  proceedings  and  journals. 


Ahbou^  tohitioa  tedmiqnet  for  both  conyressible  and 
faKonqrressible  flows  have  eviflvc^  tedmiques  for  anaNzing  these 
two  dass  of  flows  are  necessarily  different  Con^ressmle  viscous 
flows  are  governed  by  a  set  of  pmabolic  pahial  differential 
equations.  At  first  pobted  out  by  Oocco  [Re£  41  mrching  in 
time  from  an  arbitrary  set  of  initial  conditiom  is  indeed  the  most 
effidem  ww  of  solviiu  unsteady  conqirestible  viscous  flows. 
IncomprestiDle  viscous  flows  on  the  other  hand  are  governed  by  a 
mixed  set  of  elliptic  and  paisholic  differential  equations,  making 
noniterative  time  mardiing  qiproadies  difficult,  except  in  qiedal 


De^te  the  differem  mathematical  duuacteristio  of 
eomprestible  and  inomimrestfble  flows,  enough  commonalties 
exist  between  diese  two  dass  of  flows  to  encourage  adaptation  (rf 
algorithms  developed  for  incompmsible  flows  to  compressible 
flows  and  vice  versa.  For  ^  pardKdic  part  xd  the 

Incompressible  flow  equations  (eig.  the  vortid^  tranqrort 
equatira)  may  be  hitegrated  n^  sdiemes  sudi  as  the 
MacCormadc  sdieme,  or  inqilidt  approximate  factmization 
idiemes  wfaidi  are  the  traditional  tools  in  eonqrressible  flow 
solvers.  Grid  generation  sdiemes,  sod  graphia  pr^  and  post¬ 
processing  toou  are  all  quite  similar,  iriietber  one  solves  the 
cotqpffliible_flow  aquations  or  m  incompreis^  flow 


Nonetheless,  compressible  flow  solvers  have  evolved  to 
a  point  where  they  may  be  routinely  used  to  stuify  a  variety  of 
ffoblems  ranging  from  l-D  subsonic  flow  to  unsteady. 


I  fTTrTTTTTtTf-tvTl^  r 
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pressure  derivative  drops  out  in  the  steady  sute,  and  physically 
correct  solutions  are  achieved.  If  the  aim  is  to  achieve  time* 
accurate  calculadons,  either  the  artificial  pressure  derivative 


»I*Tn  r*  ■ » . 0  •Y4  I  ■  43  I  u-WXeTvpi  fTeT|TW^;^f7jVTM  n 


I  portion  of  the  inconmressible  Dow  equations  to  be  one  of  toe 
mSculties.  The  tramtional,  iterative  solution  of  the  elliptic 
portion  of  the  governing  equations  (where  the  unknown  may  be 

Ipcessure,  or  the  scalar  components  of  a  vector  stream  function) 
flooverfes  very  slowly,  and  the  convergence  rate  usually 
det^orates  at  nigh  Reynolds  numbers. 


it^  and  forces  very  small  time  steps) 
within  each  time  step  should  be  used. 


or  an  inner  iterative  loop 


us^y 


I  One  of  the  commonly  used  approadies  for  solving  2*D 
incompressible  flows  is  the  vorticiw*veloci^  or  the  vorticity* 
stream  function  formulation.  Mehta  [Ref.  S],  Wu  [Rel  6],  Thames 
[Ref.  7]  and  others  have  used  this  approach  In  Mehta’s  work, 
nnstearfy  viscous  flow  past  stationary  and  oscillating  airfoils  were 

IAudied.  He  used  an  approximate  factoriation  scheme  to 
integrate  the  parabolic  portion  of  the  governing  equations,  and  a 
direct  tet  Poisson  solver  to  invert  the  elliptic  equation  for  the 
stream  function.  He  experienced  difficulties  in  obtaining  accurate 

(solutions  at  Reynolds  numbers  as  low  as  5,000.  Some  researchers 
Wu)  have  circumvented  the  need  to  iteratively  solve  the 
Irown  equation  for  velod^  or  stream  function,  by  recasting  the 
differentia  equation  as  an  integral  relation.  The  operations  count 
tot  ♦*'<«  approach  dramatically  increases  with  Reynolds  number. 

(This  approach  requires  costly,  numerical  evaluation  of  volume 
hxtegrw  in  three-cimensions.  Thompson  [Rel  8]  studied  the  3*D 
jet-in*cross  wind  problem.  As  in  the  previous  cases,  the  analyses 
were  restricted  to  relatively  low  Reynolds  numbers. 

I  A  different  approach  for  solving  the  incompressible 

Navier-Stokes  equations  is  to  soNe  them  in  the  primitive  variables 
tana  (pressure-velocity).  Again,  a  variety  of  approaches  for 

I  solving  these  equations  are  possible.  In  some  snidies  (e.g.  Coda, 
Rel  9],  the  u*,  v*  and  w*  momentum  equations  are  integrated  in 
time  using  an  implicit  or  an  explicit  time  marching  scheme.  The 
pressure  ueld  which  appears  in  these  equations  is  evaluated  at 
eadi  time  step  by  solving  a  Poisson's  equation  for  pressure, 

I  usually  using  an  iterative  scheme.  By  using  a  suggered  storage 
sdieme  where  the  pressure  field  is  stored  at  cell  centers,  while  the 
velocity  fields  are  stored  at  face  centers,  it  is  possible  to  avoid  use 
of  nonphysical,  extraneous  boundan  condiuons  for  pressure  at 

(the  solid  surface.  Brandt  [Ref.  10],  Ta'asan  [Ref.  11]  nave  shown 
that  it  is  possible  to  accelerate  the  iterative  solution  of  the 
pressure  Geld  using  classical  multigrid  techm'ques. 

An  alternative  to  solving  a  Poisson  equation  for  pressure  is 

I  the  Marker  and  Cell  (MAQ  lugorithm  developed  at  Los  Alamos 
^  Weldi  et  aL  [Ref.  12].  At  each  time  step,  the  velocity  Geld  is 
mst  updated  using  an  explicit  time  marching  scheme,  using  the 
pressure  field  from  the  previous  time  step.  Next,  the  correction  to 

Itte  pressure  at  the  center  of  each  cell  is  evaluated  by  summing  up 
the  «"»««  flux  through  the  six  faces  of  the  cell.  If  a  ceD  is 
■rmiwiiiatitig  maw,  then  the  pressure  value  is  inaeased  to  repel 
fluid  away  from  the  celL  If  a  cell  is  losing  mass,  then  the  pressure 

I  value  is  lowered.  Any  change  in  the  pressure  field  produces  a 
oonesponding  diange  in  me  veloaty  field,  which  may  be 
computed  simply  by  evaluating  the  changes  in  the  pressim 


A  number  of  other  ^mroacbes  for  solving  incompressible 
viscous  flows  exist,  and  have  been  demonstrated  to  work  well  for 
a  limited  number  of  test  cases.  A  direct  inversion  of  the  Poisson 
equation  for  pressure  at  each  time  step  was  employed  by  Osswald 
and  Ghia  [ReL  17].  This  requires  large  amounts  memory,  because 
the  ^larse  coefficient  matrix,  that  arises  when  the  Elliptic 
emiauon  for  pressure  is  discretized,  gives  rise  to  a  full  matrix 
when  inverted.  Thus,  even  a  veiy  sparse  grid  containing  1,000 
nodes  will  require  one  million  words  of  storage  to  store  the 
inverted  matrix.  Furthermore,  a  direa  solution  to  the  eUiptic 
equation  alone  does  not  necessarily  guarantee  stable,  accurate 
solutions  of  the  entire  equation  set  Hafez  et  al.  [Ref  18]  have 
looked  at  the  decomposiuon  of  the  velocity  field  into  a  routional 
part  and  an  irrotational  part.  This  approach  shares  many  of  the 
drawbacks  of  the  tradiuonal  veloaty-vorticity  form  previously 
discussed.  Higher  order  projection  methods  for  solving  in¬ 
compressible  viscous  flows  nave  been  documented  by  Bell  et.  aL 
^et.  applied  to  low  Reynolds  number  flows 


coordinates. 

Bufldiag  upon  the  earlier  works  Chorin  [Ref  14]  and 
and  Kuuer  [Ref  IS],  Kw^  et  al.  [Ref  16]  have  developed 


flows.  In  their  scheme,  known  as  the  me - -  , - 

compressibility,  an  artificial  pressure  derivative  with  respect  to 
time  is  append^  to  the  continuiQr  equation.  The  entire  system  of 
equations  is  solved  by  a  time  marching  scheme,  as  in  a  compress- 
im  flow.  If  only  a  steady  state  is  of  intere^  then  the  added 


The  methods  for  solving  incompressible  viscous  flows 
discussed  above  have  several  drawbacks: 

a)  Most  of  them  are  only  second  order  accurate  in  space,  and 
fint  or  second  order  accurate  in  time.  Before  these  schemes  can 
be  applied  to  phenomena  such  as  direct  numerical  simulation  of 
turbulence,  it  will  be  necessary  to  raise  the  spatial  and  temporal 
accuracy  to  fourth  or  higher  order. 

b)  The  iterative  convenence  of  the  elliptic  portion  of  the 
solvers  deteriorates  at  high  Reynolds  numbers. 

c)  In  some  instances  (e.g.  method  of  pseudo-compressibility) 
a  trade  off  exists  between  temporal  accuracy  and  convergence 
speed. 

d)  These  methods  do  not  take  advantage  of  the  vast  progress 
that  has  been  achieved  in  the  solution  of  steady  viscous  Qows.  For 
example,  with  rare  exceptions,  multigrid  acceleration  of  the 
Poisson  solvers  has  not  oeen  anempted.  Acceleration  of  the 
iterative  solution  of  the  pressure  field  to  convergence  using 
spatially  varying  time  steps  and  grid  sequencing  have  also  not 
been  extensively  used.  Although  these  techniques  are  primarily 
intended  for  steady  state  solutions  in  compressible  flows,  there  is 
no  reason  why  these  strategies  can  not  be  used  to  solve  the  elliptic 
partial  differential  equation  governing  the  pressure  (or  a  vector 
stream  function). 

e)  There  has  been  a  growing  interest  in  the  use  of  massively 
|»rallel  computer  architectures  such  as  the  Connection  Machine 
to  solve  unsteady  viscous  flows.  Main  of  the  conmressible  flow 
algorithms  have  already  been  adapted  for  use  on  these  machines 
[Kef  20-22].  There  is  a  need  to  develop  new  procedures  and 
modify  existing  algorithms  for  incompressible  flows,  on  parallel 

MATHEMATICAL  AND  NUMERICAL  FORMULATION 


The  ( 
efifident,  and 


objective  of  the  present  researdi  is  to  develop 
1  accurate  solution  techniques  for  the  analysis  of  3-D, 


a)  The  sdiemes  should  be  fourth  or  higher  order  accurate 
in  Qiaee,  and  second  or  higher  order  accurate  m  time. 

b)  The  solution  techniques  should  be  capable  of  handling 


cnoplei  iotenial  and  external  flows.  That  is,  the  actions  and 
the  solution  procedures  should  be  cast  in  a  curvilinear,  tiine* 
deforming  coordinate  ^tem. 

c)  Ibe  solution  procedures  should  work  for  a  wide  range 
of  Reynolds  numbers,  with  no  appreciable  loss  in  solution 
effidenqr. 

d)  The  solution  procedures  should  be  tailored  for  efficient 

on  the  current  generation  of  vector  and  massively 
parallel  computer  architectures. 


An  coupled  system  ot  equations  for  these  delu  variables 
may  now  be  written.  For  example,  consider  the  u*  momentum 
equation  (with  densi^  assumed  to  be  unity): 

Ut  +  (U?),  +  (UV)y  +  ap/8X  -  V  (U„  +  Uyy) 

For  the  sake  of  illustration,  let  us  assume  that  a  second 
order  accuracy  in  time  is  acceptable.  Iben,  the  time  derivative 
au/at  will  be  approximated  as 

au/at  ■  (u****-u“)/At 


With  the  above  goals  in  mind,  a  solution  procedure  for 
solving  3>D  unsteady  incompressible  flows  has  been  developed. 
Jbe  key  foatures  of  tne  present  scheme  are  listed  below. 

n)  Ibe  primitive  variables  (p,u,v,w)  are  the  primary 
unknowns  in  the  present  formulatioa  Depending  on  the 
turbulence  models  to  be  use^  additional  unknowns  such  as  the 
turbulence  kinetic  energy  k,  dimipation  rate  e,  the  Reynolds  stress 
components  uV  may  need  to  be  evaluated,  fo  three-dimensional 
flows,  it  is  believed  that  solving  for  the  primitive  variables 
(p,u,v,wl  will  be  more  convenient  than  use  of  vorticity-vector 
stream  ninction. 

b)  The  present  scheme  is  iterative  in  nature.  That  is,  at  each 
time  step,  the  flow  properties  are  updated  in  an  iterative  fashion. 
Such  an  iterative  procedure  is  necessary,  because  one  of  the 
unknowns,  luunely  the  pressure  is  governed  1^  an  elliptic  PDE.  In 
some  approaches,  such  as  the  pseudo-compressibliw  approach,  a 
noniterative  marching  scheme  has  been  used.  However,  such 
schemes  trade  temporal  accuracy  for  the  ability  to  achieve  a 
convergent  steady  state  solution. 

c)  The  parabolic  portion  of  the  governing  equations  e^lqy  a 
mgh  order  (second  or  fourth  order  in  time)  implicit  time 
marehing  scheme.  Since  there  is  at  least  one  flow  variable 
(pressure)  t^t  must  be  iteratively  solved  for,  there  is  no  reason 
not  to  ufMate  the  rest  of  the  flow  variables  (u,v,w)  during  each 
pressure  iteration.  The  use  of  implicit  schemes  removes  the 
necessity  to  dioose  small  time  steps  required  by  stability 
considerations. 


The  other  terms  in  the  above  equation  will  be  evaliuted  at 
the  'll*  1/2'  time  level: 

.  (ua+U  +  oa)/2 
m  (ua+yt-i  +  n*)/2 

The  spatial  discretizations  may  be  carried  out  using  either 
a  second  order  accurate  central/upwind  difference  form  or  a 
higher  order  form.  The  higher  order  spatial  accuracy  may  be 
acnieved  on  uniform  grids  using  Fade'  ^)proximations  to  the 
derivatives;  on  highW  stretched  mds,  higher  order  accuracy  may 
be  adiieved  using  a  Lajgraneean  nt  to  the  variables  u^ ,  uv ,  p  etc 
In  high  Reynolds  number  flows,  the  Lagrangean  fit  need  not  be 
equal^  weighted  about  the  node,  but  may  be  biased  in  the 
direction  of  the  flow.  For  example,  when- the  flow  is  from  left  to 
r^t,  if  the  Lagrangean  interpolation  of  u^  is  done  using  nodes 
omy  to  the  left  of,  and  including,  the  cunent  node  then  an  upwind 
formulation  results.  The  details  of  how  the  spatial  discretization  is 
done  do  not  change  the  iterative  solution  process  being  described. 

If  the  quantities  such  as  u^  ,  uv  and  p  appearing  in  the 
above  dimerization  are  linearized  about  known  information  u° 
and  then  a  difference  equation  linking  au  ,  av  and  ap 

results.  Sudi  an  equation  is  given  tor  the  u-  momentum  equation 
below: 


au/^t-f|j^u**^‘*au)-K^u**^*^av-*-v*'*^^Ak-iau)] 


Present  Sdieme: 

For  the  sake  of  convenience,  the  details  of  the  present 
idieme  are  described  in  a  Cartesian  coordinate  system,  and  for  2* 
D  flows.  Since  the  governing  equations  may  be  cast  in  a 
curvilinear,  non-orthogonaL  time*deforming  coordinate  system  in 
a  form  ve^  rimiiar  to  the  Cartesian  form,  application  of  the 
presem  al^rithm  on  a  curvflinear  grid  is  straightforward. 

The  goal  of  the  present  scheme  is  to  advance  the  flow 
properties  (p,u,v)  from  almown  time  step 'n'  to  the  next  time  step 
v-fl'.  Let  X  be  an  iteration  counter.  Then  a  quantity  such  as 
denotes  the  variable  u  at  the  time  level  'n-t- 1'  and  iteration 
level  It*.  A  good  starting  guess  for  the  flow  variables  at  time  level 
'll*  r  at  the  start  of  the  iteration  process  is  these  variables  at  the 
previous  time  level  That  is, 

-n* 

-v« 

-P^ 


•  -  [(u*'*^^*^  •  u*)/at  *  (u*)  </uv)  * 


Here  <x  >  for  suitable,  high  order  upwind 

or  central  approximauons  to  the  spatial  derivatives. 

Note  that  the  right  side  of  the  above  equation  is  simply  the 
Crank'Nicholson  approximation  to  the  u>  momentum  equation.  If 
the  right  side  is  dnnn  to  zero,  then  the  unsteady  u-  momentum 
equation  will  be  fiilly  satisfied  at  the  cunent  time  level  n-f  1. 

A  similar  equation  inxy  be  written  for  the  v-  momentum 
equation,  linking  the  quantities  au  ,  av  and  ap.  In  the  case  of 
continuity  equaoon,  one  can  draw  upon  the  Marker  and  Cell 

airqadi,  to  link  the  iterative  changes  in  pressure  to  changes  in 
odty,  and  write 

pap- -(8,11+8  yv)***^*» 


We  also  define  "deltt  quantities'au  ,av  andap  sudt  that 

Tima,  foe  goal  of  foe  iterative  process  at  each  time  step  is 
to  drive  these  deltt  quantitiesau  ,av  smdap  to  zero. 


Herep  it  a  free  parameter,  that  may  even  vary  from  node 
to  node. 

It  should  be  noted  that  the  addition  ofp  ap  to  the  left  side 
the  above  equation  is  not  equivalent  to  a  pseudo* 
eompressibili^  qmroach.  As  long  as  ap  is  driven  to  zero,  the 
discretized  form  m  foe  continuity  equauon  is  exactly  utisfied  at 
eadi  time  step. 

.^iplying  the  above  ditcretizations  in  time  and  space  at  all 


■  tbe  nodes  in  the  flow  field,  a  ^tem  of  simultaneous'equations 
tesults  for  the  quantity  eq  equal  to  (su ,  av.  tf).  This  system  may 

Ibefonn^ywrittenas: 

lAl{aq}-{R} 

Here,  the  right  hand  side  is  the  governing  equation,  with 

I  file  tempor^  and  spatial  derivative  approximated  as  discussed 
above.  The  right  side  also  contains  the  time  derivatives  that 
qipear  in  the  governing  equation.  In  traditional  iterative  schemes 
sudi  as  the  pseudocompressibiUty  scheme,  the  right  side  contains 

■  only  the  spatial  derivatives.  Thus,  in  these  scoemes,  onlv  the 
steady  state  solution  is  guaranteed.  In  the  present  approach,  the 
time  accurate  solution  at  each  time  step  is  guaranteed,  if  the  right 
side  can  be  driven  to  zero. 

.  PRELIMINARY  RESULTS 

I  The  above  procedure  for  solving  2-D  and  3<D 

biconqiressible  viscous  flows  is  being  tested  on  several  2>D 
problems.  A  preliminary  version  of  the  procedure  for  computing 

I  external  flows  over  tlnee<dimensionaI  geometries  also  exists. 
Figure  1  shows  the  computed  surface  pressure  distribution  over  a 
body  of  revolution  tested  in  Germany,  for  which  experimentally 
measured  dau  is  arailable.  It  is  clear  that  the  propcMd  iterative 

Isdieme  gives  numerical  results  that  are  in  good  agreement  with 
csperiments. 

The  full  p^ier  wfll  give  the  following  additional  results: 

I  a)  Incompressible  flow  past  an  airfoil  at  Reynolds  numbers 

SOOOtoiiMilhW 

b)  Steady  3*D  separated  flow  past  an  ellipsoid  of 
revolution  at  an  an^e  of  attack. 

■  These  cases  have  been  chosen  because  of  the  availability 

of  good  quali^  experimental  dau ,  and  tbe  simplicity  of  the  grid 
feneration. 
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CHAPTER  I 


INTRODUCTION 


The  accurate  computation  of  three-dimensional  unsteady  incompressible  flow 
problem  is  one  of  great  interest  to  researchers  working  in  Helds  of  aerodynamics, 
hydrodynamics  and  biofluid  mechanics.  The  flow  over  complex  submarine  shapes,  flow 
past  underwater  propeller,  flow  within  turbomachinery,  and  flow  in  blood  vessels  with 
compliant  walls  are  examples  of  such  flows.  Accurate  and  efficient  computation  of  such 
flows  at  high  Reynolds  numbers  is  presently  not  possible  due  to  the  mixed  (elliptic- 
parabolic)  nature  of  the  governing  equations.  Indeed,  methods  for  three-dimensional 
incompressible  flows  lag  behind  three-dimensional  compressible  flows  by  several  years. 
Until  accurate  and  efficient  methods  for  three-dimensional  incompressible,  unsteady  flows 
become  available,  it  will  not  be  possible  to  anempt  challenging  problems  such  as  the  first 
principles  based  on  direct  simulation  or  large  eddy  simulation  of  turbulent  flows  over 
complex  geometries.  The  lack  of  such  tools  is  one  of  the  principal  reasons  that  the  first 
principles  based  prediction  of  turbulent  flows  past  and  through  complex  configurations  has 
not  been  extensively  attempted  to  date. 

As  Gresho  and  Sani  (ref.l)  pointed  out,  the  pressure  is  a  somewhat  mysterious 
quantity  in  incompressible  flows.  It  is  not  a  thermodynamic  variable  since  there  is  no 
'equation  of  state'  for  an  incompressible  fluid.  It  is  in  one  sense  a  mathematical  anefact  -  a 
Lagrange  multiplier  that  constrains  the  velocity  field  to  remain  divergence-free  ;  i.e 
incompressible  -  yet  its  gradient  is  a  relevant  physical  quantity ;  a  force  per  unit  volume.  It 
propagates  at  infinite  speed  in  order  to  keep  the  flow  always  and  everywhere 
incompressible ;  le.  it  is  always  in  equilibrium  with  a  time-vaiying  divergence-free  velocity 
Held. 


One  might  have  the  idea  that  the  compressible  Navier-Stokes  equation  solvers  can 
compute  incompressible  Hows  using  compressible  How  methods,  and  setting  the  Mach 


number  to  be  very  low.  But  this  idea  becomes  impractical  at  very  low  Mach  numbers 
because  the  compressible  Navier-Stokes  equation  solvers  have  a  singular  behavior  as  the 
Mach  number  approaches  zero.  This  leads  to  an  ill-conditioned  stiff  system  of  equations 
and  consequently  very  slow  convergence,  or  even  divergence  of  the  solution  with  time. 
.This  stiffness  can  be  explained  as  a  time  step  limitation  (ref.2).  We  note  that  all  explicit 
-methods  for  solving  the  compressible  Navier-Stokes  equations  are  limited  to  a  time  step  ' 
'which  is  less  than  that  given  by  the  CFL  condition.  For  example,  in  two-dimensions  : 

_ 1 _ 

(|ul/Ax)-t-(|vl/Ay)-i-a[(l/Ax)^  +  (l/ Ayf] 

where  a  is  the  speed  of  sound.  From  this  condition,  we  observe  that  At  approaches  zero 
as  the  speed  of  sound  approaches  infinity.  As  a  result,  an  "infinite"  amount  of  computer 
time  would  be  required  to  compute  a  truly  incompressible  flow  in  this  manner.  Implicit 
methods  will  permit  a  larger  At ,  but  the  maximum  value  is  nonnally  less  than  100  times 
that  given  by  Eq.(l.l)  because  of  truncation  errors,  approximate  factorization  errors,  and 
so  on.  Thus,  even  if  an  implicit  scheme  is  used,  it  is  not  practical  to  compute  a  truly 
incompressible  Navier-Stokes  solution  using  compressible  flow  methods. 

The  significant  difficulty  in  solving  incompressible  Navier-Stokes  equations  is  that 
the  governing  equations  are  a  mixed  elliptic-parabolic  type  of  partial  differential  equationsr 
The  continuity  equation  does  not  have  a  time  derivative  term  and  is  given  in  the  form  of  a 
divergence-free  constraint  This  is  another  major  difference  between  the  incompressible  and 
compressible  Navier-Stokes  equations.  The  absence  of  a  time  derivative  term  in  the 
continuity  equation  prohibits  time  integration  of  continuity  equation  by  a  time  marching 
scheme.  The  compressible  Navier-Stokes  equations,  on  the  other  hand, are  efficiently 
integrated  by  time  marching  schemes  because  they  are  a  set  of  parabolic  partial  differential 
equations. 

One  of  the  commonly  used  approaches  for  solving  two-dimensional  incompressible 
flow  is  the  vorticity-velocity  or  vorticity-stream  function  formulation  (ref.  3,4,5).  This  is 
very  efficient  for  two-dimensional  problems,  but  this  approach  can  not  be  extended 
straightforwardly  to  three  dimensions.  Consequently,  the  incompressible  Navier-Stokes 
equations  for  three-dimensional  problem  are  normally  solved  in  their  primitive  variable 
form  (p,u,v,w).  Most  methods  using  primitive  variables  may  be  classified  into  three 
groups.  The  first  approach  is  the  pressure  Poisson  method  or  Marker-and-Ccll  (MAC) 


method  which  was  first  introduced  by  Harlow  and  Welch  (ref.6).  In  the  pressure  Poisson 
method,  the  velocity  field  is  advanced  in  time  by  solving  the  momentum  equations  with  a 
stable  explicit  or  implicit  time  marching  scheme.  Then  the  pressure  field  is  evaluated  at  each 
time  step  by  solving  a  Poisson  equation  for  pressure  directly  (ref.7)  or  iteratively 
.(ref.8,9,10).  The  continuity  equation  is  thus  satisfied  when  the  pressure  field  is  computed 
•implicitly.  This  Poisson  equation  for  pressure  is  obtained  by  taking  the  divergence  of  the 
'unsteady  momentum  equations.  The  main  idea  of  the  MAC  method  (ref.11,12),  an 
alternative  to  solving  a  pressure  Poisson  equation,  is  that  the  pressure  field  is  updated  at 
each  time  step  by  adjusting  the  pressure  by  an  amount  proportional  to  the  negative  of  the 
velocity  divergence : 

(1.2) 

Here  the  superscripts  'k'  and  'k-T  denote  the  iteration  level,  and  p  is  a  relaxation  factor. 
Usually,  a  staggered  grid  system  (ref.6)  is  used  for  the  MAC  method,  because  such  a  grid 
does  not  require  the  specification  of  pressure  on  the  boundaries  and  does  not  produce 
unphysical  oscillations  in  the  pressure  and  velocity  Helds  due  to  the  central  differencing  of 
the  pressure  gradient  term.  The  second  approach  is  a  projection  method  (or,  fractional  step 
method  )  which  was  first  introduced  by  Chorin  (ref.  13).  At  the  first  step,  an  intermediate 
velocity  is  computed  from  the  momentum  equation  without  the  pressure  gradient  term. 
Then  a  pressure  field  is  computed  which  will  make  the  velocity  field  obtained  from  the  first 
fractional  step  divergence  free.  Finally,  a  second  fractional  step  is  performed  using  the 
pressure  field  just  computed.  The  third  group  is  the  pseudocompressibility  method 
(ref.  14, 15)  which  was  also  first  introduced  by  Chorin  (ref.  16)  primarily  for  obtaining 
steady  state  solutions.  In  this  method,  an  artificial  pressure  derivative  with  respect  to  time  is 
appended  to  the  continuity  equation.  The  entire  system  of  equation  is  solved  by  a  time 
marching  scheme  developed  for  compressible  flows,  such  as  the  approximate  factorization 
scheme  (ref.  17).  If  only  a  steady  state  is  of  interest,  then  the  added  pressure  derivative 
drops  out  in  the  steady  state,  and  physically  correct  solutions  are  achieved.  If  the  aim  is  to 
achieve  time-accurate  calculations,  either  the  artificial  pressure  derivative  should  be  kept 
very  small  (which  makes  the  equations  extremely  stiff,  and  forces  very  small  time  steps)  or 
an  inner  iterative  loop  within  each  time  step  should  be  used  (ref.  18, 19).  A  concept  similar 
to  the  pseudocompressibility  method,  known  as  the  penalty  function  method  (ref.20)  is 
widely  used  in  the  finite-element  based  incompressible  flow  solvers,  which  solves  for  p  to 
satisfy : 


Xp+VV  =  0 


(1.3) 


In  this  method,  the  pressure  gradient  term  of  momentum  equation  is  eliminated  by 
substituting  Eq.(1.3)  into  the  momentum  equation,  and  then  solving  the  momentum 
.equations  with  X  0. 

*  The  methods  for  solving  incompressible  viscous  flow  discussed  above  have  several 
drawbacks : 

a)  Most  of  them  are  only  second  order  accurate  in  space,  and  first  or  second  order  accurate 
in  time.  Before  these  schemes  can  be  applied  to  phenomena  such  as  direct  numerical 
simulation  of  turbulence,  it  will  be  necessary  to  raise  the  spatial  and  temporal  accuracy  to 
fourth  or  higher  order. 

b)  The  iterative  convergence  of  the  pressure  Poisson  solvers  deteriorates  at  high  Reynolds 
numbers. 

c)  In  some  instance  (e.g.  in  the  pseudocompressibility  method),  a  trade  off  exists  between 
temporal  accuracy  and  convergence  speed. 

d)  These  methods  do  not  take  advantage  of  the  vast  progress  that  has  been  achieved  in  the 
solution  of  steady,  viscous  flows.  For  example,  with  rare  exceptions,  muhigrid 
acceleration  of  Poisson  solvers  has  not  been  attempted.  Acceleration  of  the  iterative  solution 
of  the  pressure  field  to  convergence  using  spatially  varying  time  steps  and  grid  sequencing 
have  also  not  been  extensively  used. 

e)  There  has  been  a  growing  interest  in  the  use  of  massively  parallel  computer  architectures 
such  as  the  Connection  Machine  to  solve  unsteady  viscous  flows.  Many  of  the 
compressible  flow  algorithms  have  already  been  adapted  for  use  on  these  machines.  There 
is  a  need  to  develop  new  procedures  and  modify  existing  algorithms  for  incompressible 
flows,  on  parallel  machines. 

The  objective  of  this  study  is  to  develop  an  efticient  and  accurate  solution 
technique  for  the  analysis  of  two-  and  three-dimensional,  unsteady,  incompressible, 
viscous  flows.  The  key  features  of  the  present  scheme  are  listed  below  : 

a)  The  primitive  variables  (p,u,v,w)  are  the  primary  unknowns  in  the  present  formulation. 

b)  The  equations  and  the  solution  procedures  are  cast  into  a  curvilinear,  time-deforming 
coordinate  system  to  handle  complex  internal  and  external  flows. 

c)  An  iterative  time-marching  scheme  is  used. 

d)  The  present  scheme  is  semi-implicit  at  each  iteration  and  is  suitable  for  efficient 
execution  on  the  current  generation  of  vector  or  massively  parallel  computer  architectures. 


e)  The  solution  procedure  works  for  a  wide  range  of  Reynolds  numbers,  with  no 
appreciable  loss  in  solution  efficiency. 

f)  The  present  scheme  is  first  order  accurate  in  time  and  second  order  accurate  in  space,  but 
higher  order  accuracy  in  space  and  time  is  easily  achievable. 

Only  laminar  flow  is  considered  in  the  results  to  be  discussed  because  the  goal  of 
Hhis  study  is  to  develop  an  efficient  and  accurate  incompressible  Navier-Stokes  solver.  This 
'method  is  however  capable  of  handling  turbulent  flows  provided  a  suitable  turbulence 
model  is  used,  and  there  are  no  inherent  limitations  in  this  method  that  will  restrict  it  to 
laminar  flows. 
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CHAPTER  n 

MATHEMATICAL  FORMULATION 


In  this  chapter,  the  governing  equation  for  three-dimensional,  unsteady, 
incompressible,  viscous  flow  are  presented  in  toms  of  the  primitive  variables  (p,u,v,w)  in 
both  the  Cartesian  coordinate  system  and  a  curvilinear  non-orthogonal,  time  deforming 
coordinate  system. 


2.1  Governing  Equations  in  the  Physical  Domain 

The  motion  of  an  incompressible  viscous  flow  is  governed  by  the  conservation  of 
mass  and  momentum,  so  called  the  continuity  equation  and  the  Navier-Stokes  equation. 
Three-dimensional  unsteady,  incompressible,  laminar,  Navier-Stokes  equations  in  an 
inertial  Cartesian  coordinate  system  may  be  written  in  a  non-dimensional  form  as  follows : 

I  +  |-(E  -  E  .)  +  i(F  -  F.)  +  ^(G  -  O  .)  =  0  (2.1) 
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The  stress  terms  are  given  by 
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(2.3) 


In  Eq.(2.2)  and  Eq.(2.3),  u,  v  and  w  are  the  noimalized  (Zaitesian  components  of  velocity, 
p  is  the  normalized  pressure,  and  Re  is  the  Reynolds  number  defined  as  ; 


Re  = 


pVJL 


(2.4) 


where  p,  Y^,  L  and  H  are  fluid  density,  freestream  velocity,  reference  length  and 
coefficient  of  viscosity  (dynamic  viscosity),  respectively. 

The  governing  equation  (2.1)  is  a  mixed  set  of  elliptic-parabolic  panial  differential 
equations.  As  mentioned  before,  the  absence  of  a  time  derivative  in  continuity  equation  and 
the  absence  of  an  explicit  relationship  between  pressure  and  divergence-free  condition  on 
the  velocity  prohibit  time  integration  in  a  straightforward  manner  by  a  stable  time  marching 
scheme.  In  this  study,  the  continuity  equation  is  modified  to  directly  link  the  iterative 
changes  in  pressure  to  changes  in  velocity,  as  done  in  the  Marker-and-Cell  method. 


2.2  Governing  Equations  in  the  Computational  Domain 

If  the  above  equations  are  directly  used  on  a  (Cartesian  system  to  flow  past  complex 
geometries,  the  imposition  of  boundary  conditions  will  require  a  complicated  interpolation 
of  the  data  on  local  grid  lines,  since  the  computational  boundaries  of  complex  geomenles 
do  not  coincide  with  coordinate  lines.This  leads  to  a  local  loss  of  accuracy  in  the  computed 
solution  and  leads  to  a  complex  program.  To  avoid  these  difficulties,  a  transformation  from 
the  physical  domain  (Clartesian  coordinates(t,x,y,z))  to  computational  domain  (generalized 
curvilinear  coordinates(t,^,Ti,^))  is  used.  After  transformation  from  the  physical  domain  to 
the  computational  domain,  the  governing  equations  can  be  written  as : 


where 


and 
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with  the  contravariant  velocities  U,  V  and  W : 


Us^,  +  u§»+v^y+w^j 

V*Ti,+uTl,+vTiy+wnj 

w=Ct+uC,+v;y+w; 


Here  J  is  the  Jacobian  of  transformation 
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j  a(^.T].c) 

9(x.y.z) 

'The  quantities  ,  lit  Ct  ^  presented  if  the  grid  is  in  motion  (as  in  the  case  of  flow 
past  an  oscillating  airfoil  or  a  spinning  propeller).  These  quantities  are  given  in  terms  of  the 
velocity  of  the  grid  (x^,  ,  z^)  with  reference  to  a  stadtmary  observers : 

-  yx  -  Zt  ^2 
Ct=-XtCx“  ytCy-  ZtCz 
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H 
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CHAPTER  m 

NUMERICAL  FORMULATION 


The  numerical  procedure  for  solving  the  governing  equation  is  an  iterative  time 
marching  scheme  which  attempts  to  solve  the  discretized  form  of  equations  to  a  user- 
specified  accuracy  at  any  time  step.  Details  of  the  iterative  process  are  given  in  this  chapter. 

3.1  Grid  Generation 

The  present  method  is  a  finite  difference  scheme  which  solves  the  discretized  form 
of  the  partial  differential  equations  at  a  set  of  discrete  points  in  the  flow  fleld.  Therefore,  a 
set  of  grid  points  within  the  domain,  including  its  boundaries,  must  be  specified  before 
solving  the  governing  equations.  Such  a  body-fitted  grid  system  may  be  generated  by 
conformal  mapping,  by  algebraic  method,  or  by  partial  differential  equation  techniques.  In 
this  study,  body-fltted  C-grid  (Fig.l)  and  H-0  grid  system  (Fig.5)  are  generated  by  an 
algebraic  method  for  two-dimensional  flow  around  NACA  0012  airfoil  and  three- 
dimensional  flow  around  the  ellipsoid  of  revolution,  respectively.  For  the  three- 
dimensional  curved  duct  problem,  a  shearedAotated  Cartesian  grid  is  used. 

3.2  Grid  Motion 

In  unsteady  state  computations,  it  is  convenient  to  use  a  moving  grid  to  account  for 
the  body  motion.  The  grid  is  attached  to  the  body  and  it  rotates  or  translate  with  the  body. 
The  grid  coordinates  can  be  advanced  explicitly  by  a  first  order  time  marching  scheme  : 


=  x"  +  x“  At 

yn4.1  _  y»  ^  yn 

jjn+l  s  Z”  +  Z“  At 


(3.1) 


However,  if  only  a  pure  rotational  motion  is  considered  (say  in  a  two-dimensional  flow 
problem),  new  coordinates  of  grid  at  any  instance  in  time  can  be  simply  obtained  by  using 
the  following  relations : 


1  1 


'x'l  [cosG  -sin 01  Fx'’ 
_sin0  COS0J  z' 


(3.2) . 


where  (x,  z)  is  the  instantaneous  x,  z  values  of  the  node  and  (x',z')  is  the  x,  z  values  of  the 
node  prior  to  rotation,  and  0  is  the  clockwise  rotation  angle.  In  such  a  case  x.^  and  z^  may 
be  found  by  analytical  differentiation  of  (3.2)  with  respect  to  time  or  from  (3.1). 

3.3  Iterative  Time  Marching  Procedure 

The  goal  of  the  present  procedure  is  to  advance  the  flow  properties  (p,u,v,w)  from 
a  known  time  level  'n'  to  the  next  time  level  'n+l*.  First  of  all,  let  us  consider  the 
momentum  equation.  Since  the  momentum  equation  is  a  parabolic  type  of  partial  differential 
equation,  it  can  be  solved  using  a  time  marching  scheme  as  follows : 


i(r'-q)  +  5,E  +8,r%8,c-*  = 


+5^C» 


where  q  is  ^  of  Eq.(2.6)  excluding  the  first  row  element,  i.e.. 


(3.3) 


(3.4) 


Similarly,  E,  F,  G,  E„,  F„and  G„  can  be  also  defined.  For  example. 


uU  +  p^, 

vU  +  p^y 


wU-t-p^j 


(3.5) 


The  above  discretization  of  Eq(3.3)  is  Hrst  order  accurate  in  time  if 'm'  is  zero  or  one,  and 
second  order  accurate  if 'm'  is  set  to  1/2.  The  operators,  5^,6,,  and  6^  represent  second 

order  accurate  or  higher  order  accurate  spatial  differences.  The  higher  order  spatial  accuracy 
may  be  achieved  on  uniform  grids  using  Fade  approximations  to  the  derivatives;  on  highly 
'Stretched  grids,  higher  order  accuracy  may  be  achieved  using  a  Lagrangean  fit  to  the  flow 
variables.  In  high  Reynolds  number  flows,  the  Lagrangean  fit  need  not  be  equally  weighted 
‘about  the  node,  but  may  be  biased  in  the  direction  of  flow.  For  example,  when  the  flow  is 
from  lef.  to  right,  if  the  Lagrangean  interpolation  of  flow  variables  is  done  using  nodes 
only  to  the  left  of,  and  including,  the  current  node,  then  an  upwind  formulation  results. 

If  the  Newton  iteration  method  is  applied  to  solve  this  unsteady  flow  problem, 
Eq.(3.3)  is  rewritten  as  follows  : 

JLfqn+i.  k+1  -  q"  u  Sri"'""'-  +  5r  = 

^  f  ^  ^  ^  (3.6) 

k^l  ^  k+1  ^  g^Qn+l.  k-Kl 


Following  a  local  linearization  of  E,  F,  G,  Ey.  F  and  Gy  about  the  'n+m'  time  level  and 
at  the  'k'  iteration  level,  one  may  have 


Aq  =  ©  k 


(3.7) 


where  ©  is  a  relaxation  factor  and  A,  B  and  C  are  the  Jacobian  matrices  of  the  flux  vectors 
E  -  Ey,  F  -  Fy  and  G  -  Gy ,  respectively: 


_a(E-Ey)  a(F-Fy)  .  a(G-Gy) 

8q  ’  ■  8q  ’  aq 


and  k  is  the  residual  vector,  defined  as : 


— n+l,  k  ;rn 
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(3.9) 


Note  that  when  ^  goes  to  zero,  the  momentum  equations  in  their  discretized  form  are 
exactly  satisfied,  and  the  solution  is  independent  of  to,  and  any  approximations  made  in  the 
construction  of  A,  B  and  C. 

Next,  let’s  consider  the  continuity  equation.  As  mentioned  in  Chapter  I,  in  order  to 
.solve  incompressible  viscous  flow  problems  efficiently,  we  need  a  relationship  coupling 
xhanges  in  the  velocity  field  with  changes  in  the  pressure  field  while  satisfying  the 
'divergence-free  constraint.  In  the  present  study,  the  Marker-and-Cell  (MAC)  approach  is 
used  to  link  the  iterative  changes  between  them,  and  can  be  wiinen ; 

Ap  =  -p(VV)"'"^‘*'  (3.10) 


where  =  P  “  P 

and  ^  is  a  relaxation  factor,that  may  even  vary  from  node  to  node  using  local  time  concept. 
Again,  when  Ap  goes  to  zero,  the  continuity  equation  is  exactly  satisfied  at  each  time  step, 
even  in  unsteady  flows. 

In  curvilinear  coordinate  system,  Eq.(3.10)  can  be  rewritten  as : 


(3.11) 


The  contravariant  velocities,  U,  V  and  W  are  already  defined  in  Eq.(2.8). 

£q.(3.10)  states  that  if  a  cell  is  accumulating  mass,  then  the  pressure  value  at  next 
iteration  is  increased  to  repel  fluid  away  from  the  cell.  If  a  cell  is  losing  mass,  then  the 
pressure  value  is  lowered  to  draw  fluid.  Thus  the  pressure  field  is  iteratively  updated  along 
with  the  velocity  field  until  the  conservation  of  mass  is  satisfied. 

Combining  the  momentum  equation,  Eq.(3.7)  and  the  continuity  equation, 
Eq.(3.1 1) ,  and  applying  the  numerical  discretization  in  time  and  space  at  all  nodes  in  the 
flow  field,  a  system  of  simultaneous  equation  results  for  the  quantity  A^  equal  to 


.  This  system  may  be  fOTmally  written  as : 


[M]{A^}  =  {R} 


(3.12) 


Here,  since  the  right  hand  side  is  the  discretized  form  of  the  unsteady  governing 
equations,  as  long  as  {A^}  is  driven  to  zero,  the  discretized  form  of  unsteady  Navier- 
Stokes  equations  are  exactly  satisfied  at  physical  time  level  'n+l'. 

Although  the  matrix  [M]  is  a  sparse,  banded  matrix,  direct  inversion  of  this  matrix 
•requires  a  huge  number  of  arithmetic  operations.  A  common  strategy  in  iterative  solutions 
*of  ellipdc  equations  is  to  approximate  the  matrix  [M]  by  another,  easily  invened  matrix 
’  [N] .  The  closer  the  matrix  [N]  is  to  [M] ,  the  faster  the  iterative  convergence  of  the 
solution  at  any  time  step.  In  this  study,  matrix  [N]  contains  only  the  diagonal 
contributions  of  matrix  [M] ,  and  Eq.(3.12)  becomes  an  explicit  form  which  is  easier  to  be 
tailored  for  efficient  execution  on  the  current  generation  of  vector  or  massively  parallel 
computer  architectures  than  an  implicit  form.  This  simplicity  comes  at  the  expense  of  the 
iterative  speed.  Acceleration  of  the  iterative  process  above  is  a  major  contribution  of  this 
work  to  the  state  of  the  an. 

The  spatial  derivatives  of  convective  flux  terms  are  differenced  by  using  third  order 
accurate  upwind  QUICK  (Quadratic  Upstream  Interpolation  for  Convective  Kinematics, 
ref.21)  scheme  to  reduce  unphysical  oscillations  or  false  diffusion  for  high  Reynolds 
number  flows,  and  the  spatial  derivatives  of  viscous  terms  are  differenced  using  half-point 
central  differencing.  The  spatial  derivatives  of  continuity  equation  is  differenced  with 
central  differencing  and  a  fourth  order  artificial  damping  term  is  added  to  the  continuity 
equation  to  stabilize  the  present  procedure.  The  QUICX  scheme  is  constructed  that,  instead 
of  such  a  linear  interpolation  for  the  convective  terms  as  used  in  standard  one-sided 
differencing  schemes,  a  three-point  upstream  weighted  quadratic  interpolation  is  used.  For 
example,  let's  consider  the  convective  term  in  | -direction  which  may  be  approximated  as 
follows : 
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The  curvature  terms  (CURV)  depend  on  the  direction  of  the  contravariant  velocity  U  : 
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Fig.  3.1.  Quadratic  upstream  inteipolation 

(a)  For  U  >  0 

(b)  For  U  <  0 


3.4  Initial  and  Boundary  Conditions 

The  governing  equation  (2.1)  and  (2.5)  is  a  mixed  elliptic-parabrUc  type  of  partial 
differential  equation,  and  requires  initial  conditions  to  stan  the  calculation  as  well  as 


boundary  condition  at  every  time  step.  The  parabolic  nature  of  the  flow  ensures  that  the 
flows  will  be  independent  of  initial  conditions,  after  large  number  of  time  step. 

In  the  present  work,  the  quantities  Ap,  Au,  Av  and  Aw  are  set  to  zero  at  all  solid 
and  fluid  boundaries.  The  boundary  conditions  arc  updated  after  every  interior  points 
jupdated  during  each  iteration.  Thus  the  boundary  values  as  well  as  interior  values  are 
iteratively  advanced  from  a  time  level  'n'  to  ‘n+T. 

Initial  Conditions  : 

In  the  case  of  external  flows,  we  assume  that  the  object  is  impulsively  started  from 
rest .  Therefore,  the  uniform  freestream  conditions  are  used  as  initial  conditions.  In  the  case 
of  internal  flows,  parallel  flow  solutions  (e.g.  Poiseulle  flow  in  a  square  duct)  are  used  to 
start  the  calculations. 


Farfleld  Boundary  Conditions  : 

For  external  flow  applications,  the  farfleld  boundary  is  placed  far  away  from  the 
solid  surface.  Thus,  it  is  natural  to  specify  the  freestream  values  at  the  farfleld  boundaries 
except  along  the  outflow  boundary  where  the  extrapolation  for  velocities  in  combination 
with  P  =  P-  is  used,  to  account  for  the  removal  of  vorticity  from  the  flow  domain  by 
convective  process. 


Boundary  Conditions  on  the  Solid  Surface  : 

On  the  solid  surface,  the  no  slip  condition  is  imposed  for  velocity  components.  The 
surface  pressure  distribution  is  determined  by  solving  the  normal  gradient  of  pressure  to  be 
zero; 


dn 


=  0 


(3.17) 


Some  researchers  (ref.22, 23)  obtain  the  boundary  conditions  for  pressure  from  the  normal 
component  of  momentum  equation  at  the  wall 


^  9^u„ 

8n  Re  9n^ 


(3.18) 


where  Un  is  the  normal  component  of  velocity.  In  high  Reynolds  number  flows,  the 
viscous  stress  contribution  to  the  normal  momentum  equation  can  be  neglected  at  the  wall 
and  the  grid  point  adjacent  to  the  surface  will  be  sufficiently  fine  so  that  constant  pressure 
normal  to  the  surface  can  be  assumed.  Thus  Eq.(3.17)  is  an  acceptable  boundary  condition. 


Boundary  Conditions  on  the  Cut  and  Singular  Line  : 

Since  the  C-grid  and  the  H-0  grid  which  are  used  for  two-dimensional  airfoil 
problem  and  three-dimensional  body  of  revolution  have  a  cut  and  singularlines, 
Respectively,  special  treatment  is  needed  (see  Fig.  3.2  and  3.3).  Across  the  cut  of  the  C- 
^d  system,  flow  quantities  should  be  continuous.  The  flow  quantities  on  the  cut  can  be 


•obtained  by  averaging  the  flow  properties  from  above  and  below  the  cut.  On  the  singular 
lines  that  occur  in  a  H-0  grid  system,  the  flow  quantities  are  obtained  by  extrapolating  from 


two  adjacent  integer  points  and  then  averaging  them  azimuthally  to  ensure  that  the  flow 
quantities  are  singe-valued. 


Fig.3.2  Chit  of  the  C-grid  system 


Singularline 


Fig.3.3.  Singularline  of  the  H-0  grid  system 
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3.5  Acceleration  by  Multigrid  Technique 

Since  the  matrix  [N]  (which  is  an  approximate  to  matrix  [M]  of  £q.(3.12))  is  a 
simple  diagonal-matrix,  it  leads  to  slow  convergence  of  the  pressure  and  velocity  fields  at 
every  time  step.  Use  of  such  a  simple  diagonal  matrix  simplifies  the  inversion,  and  makes 
the  flow  solver  100%  vectorizable  and  parallelizable.  To  accelerate  the  present  procedure,  a 
multigrid  technique  (C^^arse  Grid  Correction  method)  is  applied  in  this  study. 

The  principles  behind  the  present  muldgrid  technique  are  as  follows.  The  quantities 
(Au,  Av,  Aw,  Ap)  may  be  viewed  as  Fourier  series-like  sums  made  of  components  of 
different  wave  lengths.  An  extremely  coarse  grid  linking  a  point  to  a  node  several  units 
away  is  effective  in  computing  the  long  wave  length  components.  A  very  fine  grid  is 
effective  in  computing  the  short  wave  length  components,  and  is  very  inefHcient  for 
computing  the  long  wave  length  components.  The  muldgrid  technique  attempts  to  compute 
these  individual  components  of  Aq  on  grids  of  several  levels  efficiendy.  When  the  process 
converges,  of  course,  the  discretized  equations  (i.e.  RHS  of  Eq.(3.7)  and  (3.11))  are 
exactly  sadsfied  on  the  finest  grid. 

The  coarse  grid  correction  algorithm  presently  used  (given  here  for  2-grid  sequence 
for  simplicity)  is  as  follows : 

i)  Compute  the  residual  {R}  appearing  on  the  right  hand  side  of  Eq.(3.12)  on  the  fine  grid 

_n+l.k 

using  q 

ii)  Transfer  the  residual  firom  the  fine  grid  to  the  coarse  grid  using  the  injection  operation, 
1  h  ^  •  An  injection  operation  is  given  at  any  node  (i  j)  m  two-dimensional  case  by 


+  T(Ri+l.j-l  +  Ri-l.j+l  +  Ri+l,j+l ) 


and  in  three-dimensional  case : 

”  ^ij,k  "*■  ^i+l.j.k  "*■  ^i,j+l,k 

Ri,j.k-1  +  Ri.j.k+l)  +  ^{Ri-l.j-l.k  +  Ri+l.j-l.k 

*^i+l,j+l.k  +  *^i-l.j+l.k  +  Ri.j-l,k-l  Ry+Uk-l 
*^i.j+l.k+l  +  Ri.j-l.k+1  Ri-l,j,k-l  +  Ri+l.j.k-1 


(3.19) 


(3.20) 
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+  Ri+l.j.k+1  +  Ri-l.j.k-!-l)+ g  (Ri-l.j-l.k-1  +  Ri+l.j-l.k-1 
+  Ri+l.j+l.k-1  +  *^i-l,j+kk-l  +  Ri-l.j-l.k+1  +  Ri+l.j-l.k+1 
+  Ri-l.j+l.k+1  +  ^i+l,  j+l.k+1 ) 

iii)  Compute  the  quantity  M  at  every  point  on  the  coarse  grid  by  solving  the  system  of 
equation : 

[NJlAq/Jj^llS^R}  (3.21) 

iv)  Interpolate  the  Aq  values  computed  in  step  (iii)  back  on  to  the  fine  grid  by  using  the 
bilinear  interpolation. 

v)  Compute  the  updated  values  of  the  flow  properties  as  +  Aq . 

Repeat  step  (i)  -  (v)  till  Aq  is  driven  to  zero. 

The  present  2-D  solver  accepts  grids  upm  3  levels. 

To  the  writer's  knowledge,  the  multigrid  technique  in  unsteady  incompressible 
flows  has  been  applied  only  to  pressure-Poisson  equation.  The  u-,  v-  and  w-  momentum 
equations  are  usually  solved  only  on  a  single  grid.  The  present  work  fully  exploits  the 
benefits  of  the  multigrid  method  for  all  the  equations,  while  keeping  the  form  of  the  matrix 
[N]  extremely  simple.  This  allows  use  of  larger  time  steps  and  improved  convergence  as 
discussed  on  Chapter  IV.  The  present  investigatOT  applied  a  conjugate  gradient  like  scheme, 
called  the  GMRES  (Generalized  Minimal  Residuals)  to  solve  Eq.(3.12).  The  matrix  [N] 
was  used  as  the  preconditioner.  The  success  of  the  GMRES  scheme  crucially  depends  on 

the  closeness  of  [N]  to  [Ml .  That  is  the  eigenvalues  of  the  matrix  [l  -  N”*m]  must  be 
small  and  closely  packed.  The  use  of  GMRES  with  [N]  as  a  preconditioner  was  not 
successful. 


C 


CHAPTER  IV 


RESULTS  AND  DISCUSSION 

In  this  chapter,  the  work  done  to  date  is  presented.  To  validate  the  present 
procedure,  three  cases  were  tested.  The  first  test  case  is  two-dimensional  unsteady  viscous 
flow  over  an  oscillating  airfoil.  The' second  is  three-dimensional  steady  flow  over  an 
ellipsoid  of  revolution.  The  third  is  the  flow  through  a  curved  duct.  Numerical  results  are 
presented  in  the  form  of  instantaneous  streamlines,  velocity  profiles,  vorticity  contours, 
surface  pressure  distribution,  and  aerodynamic  loads.  Streamlines  and  surface  pressure 
distributions  are  compared  with  flow  visualization  and  the  other  available  numerical  data . 

4.1  Dynamic  Stall  of  an  Oscillating  Airfoil 

The  computations  were  carried  out  for  a  sinusoidally  pitching  NACA  0012  airfoil, 
at  Re  =  S,000  and  k  =  0.5,  where  k  is  reduced  firequency  of  oscillation,defined  by 


where  D  is  the  radians  of  rotation  per  second  and  c  is  chord  of  airfoil.  The  physical 
interpretation  of  reduced  frequency  is  the  number  of  radians  of  oscillation  per  semi-chord 
length  of  travel.  This  case  has  been  previously  studied  by  Mehta  (ref.3)  at  NASA  Ames 
Research  Center  using  a  velocity-vorticity  formulation  and  its  flow  visualization  was 
carried  out  by  Werle  (ref.  24)  in  ONERA. 

After  the  flow  is  fully  developed  at  zero  angle  of  attack,  the  airfoil  is  allowed  to 
oscillate  in  pitch  through  an  angle  of  attack  range  firom  0  degree  to  20  degree  given  by 
a  s  10*(1  -  cos  t) .  Fig.l  shows  the  body-fitted  grid  around  the  airfoil  used  in  this  study. 
Hg.2  shows  the  instantaneous  streamlines  (actually,  called  particle  tracers  in  PLOT3D 
software),  velocity  profiles  and  vorticity  contours  at  selected  angle  of  attack.  Fig.3  shows 
the  surface  pressure  distribution.  In  genial,  the  streamline  patterns  and  surface  pressure 
distributions  are  in  very  good  agreement  with  flow  visualization  and  Mehta's  numerical 


results  except  that  the  present  procedure  predicts  a  little  earlier  generation  of  vortex  than 
Mehta's  method.  The  flow  visualizations  were  carried  out  with  air  bubble  in  the  water 
tunnel.  Here,  we  should  note  that  photographs  showing  air  bubble  trajectories  were  taken 
at  an  exposure  time  of  1/10  seconds.  Therefore,  in  unsteady  flow  the  air  bubble  trajectories 
.near  the  surface  of  airfoil  represent  neither  streamlines  nor  streaklines  because  the  pictures 
contain  many  paths  over  the  exposure  time.  On  the  orther  hand,  the  instantaneous 
‘Streamline  is  a  streamline  at  any  instant  of  time,  i.e.  we  assume  the  flowfield  is  frozen  at 
any  instant  of  time  and  draw  the  streamline.  In  other  words,  the  instantaneous  streamline  is 
equivalent  to  the  bubble  trajectories  with  an  infinitesimal  exposure  time.  Thus,  the  flow 
visualization  with  air  bubble  is  different  from  the  instantaneous  streamline,  and  should  be 
used  only  for  qualitative  comparison.  Fig.4  shows  the  lift,  drag  and  moment  hysteresis 
loops.  The  main  feature  of  dynamic  stall  which  is  significantly  different  from  static  stall  is 
due  to  the  generation  of  a  vortex  near  the  leading  edge.  This  vonex  passes  over  the  upper 
surface  of  airfoil,  creating  large  variations  in  the  aerodynamic  forces  and  moment.  From 
these  figures,  it  is  seen  that  the  growth  of  lift  during  the  upstroke  is  slow  and  gradual,  well 
past  the  static-stall  angle.  The  separation  region,  which  is  present  over  a  small  region  near 
the  trailing  edge  at  first,  moves  upstream  as  the  angle  of  attack  increases.  The  pitching 
moment  does  not  change  much  during  the  upstroke.  The  surface  pressure  distribution  at  an 
angle  of  attack  of  18.6  degree  shown  in  Fig.3  shows  another  pressure  peak  near  the  quaner 
chord.  This  indicates  the  leading-edge  vonex  is  already  generated,  and  this  can  be  identified 
in  Fig.2  (c).  As  the  leading-edge  vonex  moves  downstream,  the  chordwise  surface 
pressure  distribution  and  aerodynamic  forces  are  significantly  varied,  especially  during  the 
downstroke.  This  variation  may  depend  on  the  Reynolds  number,  airfoil  shapes  and 
reduced  frequency.  The  moment  stall,  associated  with  an  increase  of  negative  moment, 
begins  at  about  18.S  degree  in  the  downstroke. 


4.2  3-D  Steady  Flow  over  an  Ellipsoid  of  Revolution 

To  validate  the  capability  of  the  present  method  to  handle  three-dimensional  viscous 
flows,  the  present  procedure  was  tested  by  computing  the  flow  past  a  6:1  ellipsoid  of 
revolution  at  10  degree  angle  of  attack,  at  a  Reynolds  number  of  5,000.  Fig.S  shows  the 
body-fitted  grid  system.  Fig.6  shows  streamlines  over  the  body  surface.  There  is  a  limited 
amount  of  experimental  data  (ref.25, 26)  available  for  this  particular  configuration,  at  high 
Reynolds  number  (Re=7.2  x  10^  ).  Fig.7  shows  the  surface  pressure  distribution  on  the 
windward  and  leeward  sides  of  the  symmetry  plane,  along  with  the  experimental  data. 


Good  agreement  is  indent  everywhere  except  in  the  last  10%  of  the  body,  where  the 
present  laminar  simulation  inedicts  flow  separation,  and  a  flattening  out  of  the  pressure 
distribution. 


4.3  3*D  Steady  Flow  through  a  90*  Bended  Square  Duct 


To  validate  the  durability  to  handle  dnee-dimensional  internal  flow  problems,  the 
flow  within  a  square  duct  with  a  90-deg  bend  was  tested.  The  radius  of  curvature  of  the 
inner  wall  in  the  curved  section  is  1.8  times  of  dre  side  length  of  square  cross-section.  This 
particular  configuration  (Fig.8)  was  experimented  by  Humphrey  et  al.  (ref.27)  and 
numerically  calculated  by  Kwak  et  al.(ref.l9)  at  a  Reynolds  number  of  790  based  on  the 
average  inflow  velocity  aixl  hydraulic  diameter.  The  inflow  and  outflow  velocity  profile  are 
obtained  by  solving  the  equation  of  fully  developed  duct  flow  (ref.28) : 


9^u  ^  3^u  _  _ 

ft?  ^dx 


const 


(4.2) 


This  equation  is  a  standard  form  of  Poisson  equation  and  can  be  solved  by  ADI  scheme. 
Fig.9  shows  the  streamwise  velocity  profiles  compared  with  the  experimental  data  of 
Hunq)hrey  et  al.and  numerical  data  of  Kwak  et  al.  The  results  are  in  general  good 
agreement  with  experiments  and  the  other  numerical  solution.The  present  grid  system  is 
75x41x41.  Li  Fig.l0,  the  cross-sectional  veloci^  inofiles  are  plotted  at  8  »  30,  60  and  90 
deg.  The  top  side  and  bottom  side  of  cross-section  are  the  inside  wall  and  outside  wall, 
respectively,  bt  this  figure,  die  pair  of  vortices  and  the  seccmdary  flow  are  showi»^ese 
vortices  are  generated  due  to  the  pressure  difference  betweejn  the  higher  pressure  on  the 
outside  wail  and  lower  pressure  on  the  inside  waU.  Fig.ll  and  12  show  the  velocity 
magnitude  contours  and  the  vorticity  magnitude  omtours  at  the  three  selected  streamwise 
stations,  respectively.  The  plot  on  die  left  side  of  Hg.l3  is  a  sideview  of  streamwise 
velocity  profiles  at  y/yi/2  *  0.5,  which  is  midway  between  the  left  side  wall  and  the 
symmetry  plane  of  square  dua  and  the  right  side  plot  is  at  yly\i2  •  0,  which  is  on  the 
symmetry  plane.  The  inside  and  outside  wall  are  correspmiding  to  z  «  0  and  z  «  1, 
respectively.  Hg.l4  shows  streamwise  velocity  profiles  from  a  viewpoint  which  is  located 
at  upper  45°  in  the  xz-plane.  The  plot  at  z  »  0.25  is  corresponding  to  die  midway  plane 
between  the  inside  waU  and  die  plane  of  symmetry.  The  plots  at  z  «  0.5  and  0.75  are  on 


the  plane  of  symmetcy  and  the  midway  between  the  outside  wall  and  the  symmetry  plane, 
respectively.  Fig.  15  is  the  streamlines  viewing  from  the  three  different  viewpoints,  i.e., 
finnt,  oblique,  and  side  view.  We  can  see  the  vortex  pair  which  is  originating  from  about  6 
•0°.  Fig.l6  is  the  pressure  contours  in  the  curved  section  and  shows  the  higher  pressure 
the  outer  wall  due  to  die  centrifugal  forces. 


4.4  Acceleration  of  Flow  Solver  by  Multigrid  Technique 

The  multigtid  technique  was  irr^lemented  to  the  two-  and  three-dimensiond  solver. 
In  two-dimensional  case,  the  fine  grid  system  has  (81x41)  grid  points  and  the  coarse  grid 
system  has  the  half  of  the  fine  grid  points,  i.e.  (41x21)  grid  points,  and  the  coarsest  grid 
system  has  (21x1 1)  grid  points.  The  two  grid  system  consists  of  the  fine  and  coarse  grid 
system  (Fig.  4.1. (a))  and  the  three  grid  system  consists  of  all  of  them  as  shown  in 
Fig.4.1.(b).  Especially,  duce  grid  system  such  as  Hg.  4.1.(b)  is  called  V-cycle. 


G>arsestGrid 

(a) 


Fig.  4.1  Structure  of  muldgrid  cycle 

(a)  Two-giid  system 

(b)  Three-grid  system 


Hg.l7  shows  the  convergence  history  of  the  global  residual  (l^-nonn  of  RHS  of  Eq 

(3.12))  reduction  in  CPU  time  for  2-D  steady  flow  over  NACA  0012  airfoil  at  zero  angle  of 
attack.  Upto  40%  and  60%  acceleration  was  obtained  using  two-  and  three-grid  system. 


respectively.  The  CPU  time  is  based  on  25  iterations  at  each  time  step  on  an  IBM 
RISC/6000  workstation.  Fig.  18  shows  the  history  of  global  residual  of  2-D  unsteady  state 
for  sinusoidally  oscillating  airfoil  (SO  iteradons/time  step),  where  the  three-grid  system  is 
used  for  muldgrid.  The  residual  by  the  multigrid  technique  maintains  lower  level  than  that 
of  single  grid  iteration  procedure  indicating  that  the  discretized  equations  are  solved  to 
much  high  levels  of  accuracy  using  the  multigrid  technique.  The  surface  pressure 
distribution  and  dynamic  stall  hysteresis  is  nearly  the  same  as  those  of  single  grid  system 
and  are  not  plotted  here.  In  three-dimensional  case,  the  multigrid  technique  was 
implemented  to  the  flow  solver  for  90°  bended  square  duct  problem.  The  three  level  of  grid 
system  consists  of  (65x21x21),  (33x11x11),  and  (17x6x6).  Fig.l9  shows  the  comparison 
of  convergence  history  with  and  without  multigrid  method.  Here  we  got  much  better 
quality  of  solutions  with  multigrid  technique  than  were  without  multigrid  technique. 
Furthermore,  The  grid  system  (65x21x21)  is  so  coarse  that  it  can  not  detect  the  sufficiently 
strong  vortex  core  and  the  residual  of  solution  without  multigrid  technique  remains  of  the 
order  of  10°.  Thus  the  solution  with  single  grid/non-multigrid  version  is  not  accurate.  The 
comparison  of  solutions  with  multigrid  and  without  multigrid  is  shown  in  Fig.20.  The 
solutions  of  fine  grid(7 5x4 1x41)  system  with  single  grid  are  also  plotted  to  compare  the 
intensity  of  vortex.  From  this  figure,  it  is  clear  that  the  multigrid  analysis  is  adequately 
resolving  the  counter-rotating  vonex  pair. 


4.5  3-D  Steady  Flow  around  a  Marine  Propeller 

The  flow  around  a  marine  propeller  is  a  challenging  problem  because  the  geometry 
is  much  more  complex  than  the  aircraft  wing  and  helicopter  blades.  The  high  twist,  low 
aspect  ratio,  close  proximity  of  other  blades  and  high  rotational  speed  make  the  flow 
around  a  propeller  highly  three-dimensional  and  complex,  featuring  centrifugal  forces, 
formation  of  curved  tip  vortices  and  leading  edge  vortices. 

Technique  for  efficient  and  accurate  prediction  of  3-D  incompressible  viscous  flow 
around  a  marine  propeller  are  necessary  for  accurate  prediction  of  performance.  Moreover, 
the  lack  of  such  tools  is  a  major  obstacles  to  the  accurate  calculation  of  cavitation  and 
propeller  noise. 

Numerical  methods  to  solve  flow  around  propellers  range  from  Goldstein's  strip 
theory  (1929)  to  Navier-Stokes  equation  solvers.  The  Goldstein's  strip  theory  models  the 
propeller  by  a  lifting  line  vortex  in  a  potential  flow  and  assumes  that  the  wake  is  a  rigid 


helical  vorte}(  sheet  This  theory  can  handle  only  a  straight  blade  without  a  nacelle.  Sullivan 
(1977)  and  Egolf  (1979)  improved  this  theory  to  account  for  blade  sweep  and  nacelle  by 
using  the  curved  lifting  line  concept  and  vortex  ^laments.  A  review  of  potential  flow 
method  applicable  to  propellers  is  well  described  by  Kerwin  (1986).  Jou  (1982)  has 
applied  the  full  potential  equation  with  a  finite  volume  approach  to  solve  propfans  but  his 
method  can  not  catch  the  strong  rotational  flow  effects  near  the  leading  edge.  Euler 
equations  have  been  applied  by  many  researchers  such  as  Chaussee  (1979,  1983),  Bober 
(1985, 1986)  and  Whitfield  (1987).  For  a  more  accurate  prediction,  Srivastava  and  Sankar 
(1990)  developed  an  iterative  method  which  couples  the  Euler  equation  solver  and 
NASTRAN  to  model  structural  deformation  due  to  aerodynamic  forces  and  centrifugal 
forces.  Full  Navier-Stokes  equations  have  recently  been  applied  to  advanced  propfans  by 
Matsuo  (1988)  and  Hall  (1991).  Lim  and  Sankar  (1993)  extended  the  Euler  equation  solver 
of  Srivastava  and  Sankar  to  full  Navire-Stokes  equations  using  the  Roe  upwind  scheme. 
These  Navier-Stokes  equation  solvers  are  based  on  the  compressible  Navier-Stokes 
equations  and  can  not  accurately  predict  incompressible  flow  solutions.  Kim  (1989)  applied 
incompressible  Navier-Stokes  equations  in  the  cylinderical  coordinate  to  an  infinite-pitch 
marine  propeller  with  rectangular  blades  by  using  the  SIMPLER  algorithm.  Although  their 
work  can  simulate  marine-propeller  flow  fields,  including  propeller  loading  and  complex 
blade-to-blade  flow,  the  infinite  pitch  propeller  with  rectangilar  blade  shape  is  not  realistic 
and  no  experimental  data  is  availble  for  comparison. 

The  above  procedure  for  solving  3-D  unsteady  incompressible  viscous  flows 
without  cavitation  has  been  applied  to  the  flow  around  a  2-bladed  SR7L  propeller  as 
shown  in  Fig.l.  The  present  scheme  is  time  accurate  and  the  steady  state  solutions 
are  obtained  as  asymptotic  solution  of  the  time  marching  process.  Fig.2  shows  the 
H-0  grid  for  a  2'bladed  SR7L  propeller.  Fig.3  shows  the  pressure  coefficient 
distribution  at  some  selected  spanwisw  location  at  a  nondimensional  time  of  0.3  in  a 
single  grid  system  (without  accleration  by  multigrid  technique),  compared  with 
experimental  data  by  Bushnell  (1988)  and  compressible  Navier-Stokes  equation 
solutions  by  Lim  and  Sankar.  A  fairly  good  agreement  with  experimental  and  other 
numerical  data  was  achieved  except  near  the  leading  edge  region.  These 
discrepancies  are  because  at  the  nondimensional  time  =  0.3,  the  flow  has  not 
evolved  enough  to  generate  leading  edge  vortices.  At  later  time  levels,  it  is 
anticipated  that  the  suction  peak  near  the  leading  edge  will  be  higher. 
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CHAPTER  V 

CONCLUDING  REMARKS 


An  accurate  and  efficient  iterative  time  marching  procedure  for  two-  and  three-dimensional 
unsteady,  incompressible,  viscous  flow  has  been  developed.  It  has  been  applied  to  the  following 
cases  with  success : 

(a)  Massively  separated  flow  over  oscillating  airfoil, 

(b)  Three-dimensional  flow  past  an  ellipsoid  of  revolution, 

(c)  Three  dimensional  flow  through  a  square  duct  with  a  90-deg.  bend 

(d)  Three-dimensional  flow  around  a  marine  propeller. 

Good  agreement  with  published  experimental  and  numerical  data  has  been  obtained  After 
the  validation  of  the  present  procedure,  techniques  for  acceleration  were  explored.  It  was  found 
that  the  muldgrid  technique  was  efficient  in  reducing  the  CPU  time  needed  for  the  simulation  and 
improved  the  solution  quality  because  of  the  lower  residuals  achieved. 

This  report  is  a  draft  copy  of  the  Ph.  D.  dissertation  of  Mr.  Wam-Gyu  Park,  and  will  in  its 
completed  form  include  changes  suggested  by  the  thesis  committee.  A  copy  of  the  finished  form  of 
the  thesis  will  be  mailed  to  the  sponsor  around  March  1993. 


Figure  1.  Bo(fy*Fitted  Grid  Around  a  NACA  0012  airfoil 


(a)  a  =  0  (deg) 

Figure  2.  Instantaneous  Streamlines,  Velocity  Profiles,  and  Vorticity 
Contours  at  Selected  Angle  of  Attack  with  Experimental  Flow 
Visualizations 


Figure  2.  Continued. 


(c)  a  =  19.8  (deg) 
Figure  2.  Continued. 
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Figure  2.  Continued. 


Figure  2.  Continued. 
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Figure  4.  Dynamic  Stall  Hysteresis  Loops  for  a  NACA  0012  Airfoil 
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Figure  5.  Body*Fitted  Grid  Around  an  Ellipsoid  of  Revolution 


Figure  6.  Streamlines  over  the  Ellipsoid  of  Revolution 


(b)  In  The  Leeward  Plane 

Figure  7.  Surface  Pressure  Distributions  Compared  With  Experimental 


Hguie  8.  Body-Atted  grid  within  a  square  duct 
with  a  90°  bend 
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Figure  9.  Streamwise  Velocity  Profiles 


Figure  12.  Vorticity  magnitude  contours  at  three 

streamwise  stations  in  the  curved  section 
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Figure  14.  Streamwise  velocity  profiles  at  z  =  0.25,  z  =  0.5, 
and  z  =  0.75 


Figure  17.  Comparison  of  convergence  history  of 
NACA  0012  airfoil  for  steady  state 
(at  zero  angle  of  attack) 
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Figure  1 8.  Comparison  of  global  residual  history 
of  a  sinusoidally  oscillating  airfoil 
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Figure  19.  Comparison  of  convergence  history  of  3-D 
steady  flow  through  a  90°  bended  square  duct 
(Grid  :  65x21x21) 
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Figure  20.  Comparisons  of  streamwise  velocity  profile 
with  and  without  multigrid  technique 


Figure  22.  Body-fitted  H-0  grid  around  a  SR7L 
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Figure  23.  Chordwise  pressure  distribution  at  selected 
spanwise  location  ( Advance  ratio  =  0.881, 
P  =  31.0,  Re=50,000) 
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Figure  23.  continued 
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Figure  23.  continued 
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Figure  23.  continued 
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Figure  23.  continued 
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